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Appendix  A 


Previous  Truth  Model  and  Filter  Equations 
(Ref  5:9-29) 


Having  assumed  xpeak  and  ypeaJc  to  be  the  sum  of  Xp  and 
xA,  and  yD  and  yA  respectively,  the  previous  study  further 
modelled  the  dynamic  components  mathematically  as  (Ref  5:10) 


Vfc)  - 

-  1/td  XpCt)  +  wx(t) 

(A-l) 

yD(t)  * 

"  1/tD  yD(t)  +  w2(t) 

(A-2) 

where 

E{w1(t)} 

=  e{w2  (t) }  =  o  yt 

(A- 3) 

Etv^Ct)  Vj^s)}  =  E{w2(t)  w2(s)}  ■  2c^/td  6  (t-s) 

(A-4) 


E{w^(t)  w2(s) >  =0  V  t,s  (A-5) 

and  td  *  correlation  time 

w1(t),w2(t)  «  independent,  white,  Gaussian  noise 
processes 

2 

oD  *  desired  sigma  or  root  mean  square 
(RMS)  value  on  the  outputs  Xp  and  yD. 


The  atmospheric  jitter  components,  xA  and  yA  were  modelled 
as  the  outputs  of  third  order  shaping  filters  driven  by 
white,  Gaussian  noise.  This  shaping  filter  is  represented 
in  the  frequency  domain  by  the  following  transfer  function 
(Ref  5:10,12) : 


Z 

w 


Kab2 

(s+a) (s+b) 2 


(A-6) 


.1 


where  a  =  14.14  rad/sec 

b  ■  659.5  rad/sec 

K  ■  gain  adjusted  to  obtain  desired  root  mean 
square  (RMS)  jitter  output 

and  E{w (t) }  «  0  (A-7) 

E{w(t)  w (s ) }  =  1  6  (t-s)  (A-8) 

Incorporating  these  models  into  state  space  notation  yielded 
a  time  invariant  form  (Ref  5:11,13)  for  the  truth  model 


where 


and 


-  Lr^t)  +  G^t)  (a-9) 

It  =  8  by  8  truth  model  plant  matrix 

«  truth  model  state  vector  with  (8  components) 

®T  *  truth  model  input  matrix  (8  by  4) 

w,p  *  vector  of  white  Gaussian  noise  inputs  (4 
by  1) 


Efw^t)}  -  0 


(A-10) 


Efv^Ct)  wj(s)>  -  Q,p  6  (t-s) 


(A-ll) 


(A-12) 


•2 


Mercier  had  numerical  problems  with  the  standard  phase 
variable  form  for  A-9  but  successfully  used  the  Jordon 
canonical  form  as  follows  (Ref  5:13-14): 


2r<t) 


o 

G, 


0 

0 

0 

0 


0 

0 

0 

0 

1 

0 

0 

0 


0 

0 

0 

0 

0 

G, 


(A-13) 


where 


KabV(a-b) 


G2  *  "G1 


G3  -  KabV(a-b) 


and  the  output  matrix  is  then  described  by 


Sjtt) 


m 

11100000 

yp«^t(t). 

00001110 

(A-14) 
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The  state  transition  matrix  corresponding  to  A-13  is  block 
diagonal  due  to  the  independence  of  the  two  channels.  Each 
4  by  4  diagonal  block  is  equal  to  (Ref  5:14-15) 


^Txtti+l'ti) 


$Ty (ti+l,ti) 


e-At/xD 


.-aAt 


0  e"bAt  Atc“bAt 

„  «  ~-bAt 


(A-15) 


So,  an  equivalent  discrete  time  model  of  the  time  propaga¬ 
tion  of  Xq,(t)  is  given  by  (Ref  5:15) 


Sp^i+l^  *  +  — Td^i^ 


(A-16) 


where 


2Td(ti>  *  f 


^rtti+l,T)  SjtTjWg.t'OdT  (A-17) 


as  commonly  used  but  wTd  is  more  properly  defined 

/fci+l 

5rtti+l»T)  §t(t)  d£(T)(R*f  8:171) 

t .  , 


(A-17a) 


ElSra(ti>1  '  0 '  E{ara(ti,2Td(tj)}  *  9ra5ij  (A-18) 
r*1*1  t  t 

flpd  “  J  — -T  ^ fci+l * ^)  — T  (x )  — T  ( t) S.T  ( T ) ip  ( ti+i » T )  dx 

'  t j  /ml  a\ 


(A-19) 


For  simulation  implementation,  A-16  was  used  in  the  follow¬ 
ing  form  (Ref  5:16) 


2&P(ti+l} 


♦jj,  ^(t^)  +  V^Td  ^(ti) 


(A-20) 


where  is  a  Cholesky  square  root  matrix  of  QTd  such 

that 


V§^  V^ST  -  Sra  (A-2D 


and  is  a  vector  of  independent  Gaussian  noises  whose 

terms  each  have  unit  variance. 

The  truth  model  measurement  equation  expresses  mathe¬ 
matically  that  each  pixel  output  is  the  sum  of  target 
intensity  radiation,  PLIR  generated  noise,  and  background 
noise  as  follows  (Ref  5:17) : 


zjktti) 


T~  If  Imax  exp  I — ^7  f<x- 

Ap  jkth  max  (  2 Cg  L 


jkth 

pixel 


♦  I  dxdy  +  *  fjk(ti>  <A-22> 

where  ®jk^ti^  *  outPut  of  Pixel  at  time  t^ 

Ap  ■  area  of  a  pixel 

(x,y)  ■  coordinates  of  any  point  in  the  pixel 
array 

I  .  ■  maximum  target  intensity 

max 

og  -  dispersion  of  target's  Gaussian 
intensity  function 


njk(ti)  »  background  noise  for  the  jkth  pixel 

fjj^tti)  »  FLIR  noise  term  for  the  jkth  pixel 

He  assumed  that  n^Ct^)  and  were  independent  wide 

band  noises  which  could  be  appropriately  modelled  as  zero 
mean,  white  Gaussian  processes  with  the  following  statis¬ 
tics  (Ref  5:18): 

E{nkl(ti)  nkittj)>  *  al  5ij  (A-23) 

'  «t  4ij  <A-24> 

Ckl  ?  mn)  (A-25) 

Captain  Mercier  used  a  four  state  extended  Kalman 
filter  -  one  state  for  target  dynamics  representation  and 
one  state  for  atmospheric  jitter  representation  in  each  of 
two  directions  in  the  FLIR  focal  plane.  Each  state  was 
modelled  as  a  stationary,  first  order,  Gauss-Markov  process 
expressed  as  (Ref  5:19): 

xCt)  -  -  \  x (t )  +  w (t )  (A-26) 

5_2 

with  E{w(t)}  -  0  »  E{w(t)wts) }  ■  6(t-s)  (A-27) 

where  t  *  correlation  time 

w(t)  »  white  noise  process  of  strength  q 
2 

a  ■  desired  variance  of  x(t) 
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The  state  space  form  where  Xp  *  ^  XA  yD  yA^  is 

(Ref  5:20) s 

Xp(t)  «  Pj,x  (t)  +  Wp(t) 

where  0  0  0  1 


(A-28) 


-1/T, 


(A-29) 


-1/x, 


Wp(t)  «  input  white  noise  vector  of  dimension  4 


such  that 


EtWp  (t)} 


(A-30) 


E{Wp(t)Wp(s)  } 


Q<5  (t-s) 


qD  o  o  o 

0  qA°  0  6  (t-s) 

0  0  qD  0 

0  0  0  qA 


(A- 31) 


The  first  order  system  approximation  in  A-28  and  A-29  for 
the  atmospheric  jitter  described  in  A-6  is  justified  based 
on  the  discrepancy  between  the  two  break  frequencies  and  the 
dominant  effect  of  the  low  frequency  pole. 

With  conversion  to  an  equivalent  discrete  time  form  of 
equations  (A-28)  to  (A-31),  the  filter  propagation  equations 
are  (Ref  5:20-23): 


2p(ti+l)  "  ^Ftti+l,ti)  2p(ti) 


(A- 32) 


^F(ti+1) 


^F  * ti+l '  fci *  * ti  *  * ti+l '  ti * 


/*i+l  T 

+  J+  ^Flti+1'  T)^(T)^(ti+1,T) 

'  tj 


dx 


(A-33) 


where  ^(t,^) 

-(t-t^/Xp 


0 

0 

0 


0 

-(t-t^/Tj 

0 

0 


0 

-(t-t^/Xp 


0 

0 

-(t“V/TA 


(A- 34) 


and  Qp  is  a  constant  matrix  whose  entries  represent  the 

appropriate  strengths  of  the  continuous  time  noise  processes 

2 

in  (A-29)  to  obtain  the  proper  output  variances  of  oZ.  and 

2  u 
cA  (Ref  8:179).  Therefore 


2f 


!!2 

tD 

0 


ta 

0 


fs 

td 

0 


20 


and  the  second  term  in  A-33  becomes 

-2At/tp 


»D(1-* 


0 

0 

0 


0|(l-e 


■) 


0 

0 


0^(1-. 


0 

-2At/TD 


) 


(A-35) 


0j(l-e 


0 

0 

0 

-2At/x, 
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Similar  to  A-22,  the  measurement  relation  is  (Ref  5:24) 


zjk(ti) 


A  //, 
AP  ijth 
pixel 


exp 


+  (y-ypeak(ti>)2]S  ***  +  njk(ti>  +  cjk(tl>  (A-36» 


where  the  only  difference  between  A-22  and  A-36  is  that 

xpeak  and  ypeak  are  functions  of  Xp^).  Since  an  8  by  8 
tracking  window  is  used  the  measurements  are  expressed  as  a 
6 4 -dimensional  vector  and  are  assumed  to  be  a  nonlinear 
function  of  the  filter  states  plus  noise  (Ref  5:25), 


=  h^ltj.)  f1 

tL)  +  v.(ti) 

(A-37) 

where 

E{v(ti) }  = 

0 

(A-38) 

E{v(ti)vTCtj)>  = 

*F  *  6ij 

(A-39) 

Thus,  the  usual  form  of  the  extended  Kalman  filter  update 

equations  would  be  (Ref  5:26) 

Ktt^  - 

£ltJ)HT(ti)  [HCt1)P(tJ)HT(ti)  +  R(t±) 3”1 

(A-40) 

x(t^) 

*  x(t^)  +KCti){£(ti)  -  h(x(tT)  ,ti) } 

(A-41) 

P(tJ)  -  P(t^)  -  K(ti)H(ti)P(t“) 

(A-42) 

where 

3h(x,t) 

A  M 

(A-43) 

X  -  Xtt^ 

and  rtt^  is  a  64-dimensional  vector  of  actual  (based  on 
the  truth  model)  measurements .  To  avoid  the  inversion  of 
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a  64  by  64  matrix  as  required  by  A-38,  Captain  Mercier  chose 
to  implement  another  form  of  the  update  equations  known  as 
the  inverse  covariance  form  (Ref  8:238-241)  shown  below: 

P(t+)  -  l£._1CtT)  +  HT Cti> R—1  (tj^) H (ti)  ]—1  (A-44) 

E(ti)  -  P(ti)HT(ti)R"1(t.)  (A-45) 

x(tj)  -  xCt^  +  K(t±)  (£ (t±)  -  h(x(t7)  ftA)]  (A-46) 

This  set  of  equations  requires  the  inversion  of  two  4  by  4 
inatr^ce8*  Because  the  pixel  noises  are  assumed  to  be  inde¬ 
pendent,  the  Rp  is  diagonal  and  its  inverse  simplifies  to 
(Ref  5:27)  : 

Bp1  (tA)  *  (j^-)  I  (A-47) 

which  can  be  precomputed  once  offline. 

The  partial  derivative  matrix  in  equation  A-41  is 
easily  computed,  using  the  fact  that  the  integral  is  a 
linear  operator,  and  results  in  a  64  by  4  matrix  (Ref  5: 
27-28).  The  particular  components  of  this  matrix  can  be 
seen  explicitly  in  subroutine  MEASF  in  Appendix  F. 


Real  Data  Analysis  for  Spatial  Correlations 
of  Background  Noise 


Real  digital  FLIR  data  has  been  analyzed  to  establish 
the  spatial  character  of  the  background  noise.  The  data 
was  from  a  real  exercise  of  tracking  an  air-to-air  missile 
conducted  in  the  Spring  of  1978  by  the  Air  Force  Weapons 
Laboratory  (AFWL) .  The  data  from  AFWL  consisted  of  a  matrix 
of  digital  data  printed  every  half  second  of  the  test.  Each 
frame  or  matrix  of  data  is  called  a  record  with  record  num¬ 
bers  running  from  113  to  190  (38.5  seconds  in  total  dura¬ 
tion)  .  As  the  missile  is  flying  toward  the  sensor  through¬ 
out  the  exercise,  the  digital  values  toward  the  end  of  the 
time  span  represent  the  sum  of  2-by -2  or  4-by-4  arrays  of 
pixels  instead  of  individual  pixels  to  maintain  the  missile's 
image  in  the  limited  data  field  of  view.  The  background  ini¬ 
tially  is  clear  sky  with  some  ground  clutter  appearing  in 
the  latter  segment  of  frames  when  the  missile  motion  carried 
it  below  the  horizon. 

The  measure  of  correlation  between  two  sets  of  data, 

X  and  Y,  of  equal  dimension,  is  most  commonly  expressed  by 
a  correlation  coefficient,  r.  The  standard  equation  for 
the  computation  of  r  is  (Ref  11:565): 


r  « 


”  Vi  -  ( "  xi)  ( ” 
i-1  1  x  \j-l  V  \j«l 


n 

n  £ 
i-1 


n 

£ 

i-1 


(B-l) 

Using  this  relation,  the  values  in  selected  rows  and  columns 
of  the  real  digital  data  which  appeared  to  be  unassociated 
with  the  real  target  image  in  the  frame  were  evaluated  for 
their  correlation  coefficient.  For  example,  m  values  in  a 
given  row  were  evaluated  for  a  horizontal,  first  neighbor 


correlation  coefficient  by  letting  in  B-l  equal  the  first 
m-1  values  in  the  row  and  equal  the  corresponding  values 
displaced  one  pixel  to  the  right  (n*m-l) . 

Tables  B-I  through  B-IV  contain  the  correlation  eval¬ 
uations  that  were  conducted.  Table  B-I  shows  evaluations 
of  pixel  correlations  for  pixels  side  by  side  (first  neigh¬ 
bor)  in  the  horizontal  direction.  The  data  shown  is  record 
or  frame  number,  the  row  or  column,  how  many  pixels  are 
summed  to  get  the  digital  data  value,  and  the  correlation 
coefficient  (r) .  Note  that  the  rows  are  labeled  A  and  B 
to  reflect  that  the  rows  of  data  have  a  biased  average  value 
when  compared  to  the  row  above  or  below  it,  as  was  first 
discussed  in  Chapter  I.  Table  B-II  contains  correlation 
coefficients  of  each  digital  value  with  its  second  hori¬ 
zontal  neighbor.  Table  B-III  shows  vertical  first  neighbor 
correlations  evaluated  by  considering  a  column  of  data. 

When  the  data  values  represent  a  single  pixel  output  as  is 
the  case  of  the  first  six  entries  of  Table  B-III,  the  first 
neighbor  vertical  correlations  represent  an  attempt  to  com¬ 
pare  A  values  and  B  values.  Table  B-IV  represents  second 
neighbor,  vertical  correlations.  Its  first  six  entries, 
(single  pixel  values)  in  contrast  to  Table  B-III,  represent 
correlations  of  A  values  with  A  values  and  B  values  with 
B  values. 

There  are  several  notable  trends,  some  of  which  were 
unexpected,  in  this  data.  First,  horizontal  neighbors  cor¬ 
relations  using  one-pixel  values  averaged  .323  with  a 
standard  deviation  of  .119  which  can  be  well  matched  by  an 
exponential  decay  model  with  a  correlation  distance  of  one 
pixel;  similar  correlations  using  four  and  sixteen  pixel 
values  averaged  .249  (last  six  entries  of  Table  B-I)  with 
a  standard  deviation  of  .224.  Second,  horizontal  neighbors 
showed,  surprisingly,  a  slight  negative  average  correlation 
of  -.082  for  single  pixel  values  with  a  standard  deviation 
of  .113.  For  four  and  sixteen  pixel  values,  second  neighbor 
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TABLE  B-I 

First  Horizontal  Neighbor  Correlations 


No.  of  No.  of 

Entries  Pixels  Correlation 


Record 

Row 

per  Row 

per  Entry 

Coef  (r) 

114 

A1 

28 

1 

.27139 

114 

B1 

28 

1 

.30296 

114 

A12 

28 

1 

.48088 

114 

B12 

28 

1 

.35063 

119 

A1 

30 

1 

.29152 

119 

B1 

30 

1 

.43761 

119 

A12 

30 

1 

.08746 

119 

B12 

30 

1 

.36134 

181 

Al 

24 

4 

.32512 

181 

B1 

24 

4 

-.07193 

181 

B15 

24 

4 

.5637 

181 

A16 

24 

4 

.1393 

185 

Al 

24 

16 

.39189 

185 

All 

24 

16 

.1452 

One  pixel  average  .323 

standard  deviation  . 119 

Four  and  sixteen  pixel 

average  .249 

standard  deviation  .224 

horizontal  correlations  were  low,  also,  averaging  .123  with 
a  standard  deviation  of  .163.  As  mentioned  before  the  ver¬ 
tical  correlations  were  confused  by  the  A  and  B  row  biased 
averages.  In  Table  B-III  when  A  and  B  values  were  compared, 
negative  correlations  averaging  -.193  resulted  with  a 
standard  deviation  of  .166.  However,  when  A  and  B  values 
are  summed  in  the  four  and  sixteen  pixel  digital  values, 
the  average  correlation  coefficient  was  .355  with  a  stand¬ 
ard  deviation  of  .198.  Second  neighbor  vertical  correla¬ 
tions  for  single  pixel  values  in  Table  B-IV  were  very 
»  inconsistent  averaging  .160  (three  coefficients  of  twelve 
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TABLE  B-II 


Second  Horizontal  Neighbor  Correlations 


Record 

Row 

No.  of 
Entries 
per  Row 

No.  of 
Pixels 
per  Entry 

Correlation 
Coef  (r) 

114 

A1 

28 

1 

-.2266 

114 

B1 

28 

1 

-.06488 

114 

A12 

28 

1 

.03753 

114 

B12 

28 

1 

-.15656 

119 

A1 

30 

1 

-.117 

119 

B1 

30 

1 

-.1955 

119 

A12 

30 

1 

.0948 

119 

B12 

30 

1 

-.0283 

181 

A1 

24 

4 

.1475 

181 

Bl 

24 

4 

-.0197 

181 

B15 

24 

4 

.1634 

181 

A16 

24 

4 

.0934 

185 

A1 

24 

16 

.4052 

185 

All 

24 

16 

-.0513 

One  pixel  average  -.082 

standard  deviation  .113 

Four  and  sixteen  pixel 

average  .123 

standard  deviation  .163 

were  negative)  with  a  standard  deviation  of  .261.  For  four 
and  sixteen  pixel  digital  values,  the  second  neighbor  ver¬ 
tical  average  was  .335  with  a  standard  deviation  of  .155. 

Physically,  the  horizontal  correlations  were  expected 
to  be  greater  than  the  vertical  correlations  due  to  the  hori 
zontal  mechanical  scanning  process  of  the  pixels.  This 
could  be  verified  by  an  experiment  where  a  well  defined 
background  is  viewed  by  the  FLIR  in  its  normal  orientation 
and  then  viewed  by  the  FLIR  turned  90  degrees.  Spatial 
analysis  of  data  from  both  angles  would  be  able  to  sepor- 
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TABLE  B-III 


First  Vertical  Neighbor  Correlations 
A  to  B 


Record 

Column 

No.  of 
Entries 
per  Column 

No.  of 
Pixels 
per  Entry 

Correlate 
Coef  (r) 

114 

1 

24 

1 

-.51563 

114 

2 

24 

1 

-.16357 

119 

1 

24 

1 

-.1850 

119 

2 

24 

1 

-.0422 

119 

29 

24 

1 

-.0989 

119 

30 

24 

1 

-.1511 

181 

23 

31 

4 

.6059 

181 

24 

31 

4 

.13625 

185 

1 

22 

16 

.41517 

179 

21 

28 

4 

.44834 

179 

22 

28 

4 

.17058 

One  pixel  average  -.193 

standard  deviation  .166 

Four  and  sixteen  pixel 

average  .355 

standard  deviation  .198 

ate  the  contribution  to  the  horizontal  correlations  from 
the  scanning  process. 

The  correlations  for  this  data  sample  were  not  as  wide¬ 
spread  as  anticipated.  This  is  possibly  explained  by  the 
fact  that  the  background  was  clear  sky  when  the  digital 
values  were  single  pixel  outputs.  During  the  last  segment 
of  the  test  when  four  and  sixteen  pixel  sums  were  used,  the 
background  was  probably  ground  clutter  as  the  missile  had 
moved  below  the  horizon.  This  fact  may  supply  some  under¬ 
standing  of  why  spatial  correlations  for  the  last  segment 
seemed  equal  to  or  larger  than  those  of  the  single  pixel 
correlations  (i.e.  four  times  the  initial  correlation  dis¬ 
tance,  or  greater).  Obviously,  more  data  analysis  is  needed 
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Appendix  C 

Real  Data  Analysis  for  Temporal  Correlations 
or  Background  Noise 

Real  data  was  analyzed  using  REFERENCING ,  a  FORTRAN 
software  package  supplied  by  AFWL,  to  substantiate  the 
existence  of,  and  to  determine  numerical  values  for,  the 
temporal  correlations  of  the  background  noise  of  the  FLIR 
pixels.  REFERENCING,  described  in  great  detail  in  Reference 
12,  was  written  in  order  to  provide  a  computer  simulation 
capability  for  the  digital  correlation  trackers  using  real 
data  values  stored  on  magnetic  tape.  The  data  tapes  con¬ 
taining  the  raw  data  to  be  read  into  REFERENCING  were  sup¬ 
plied  by  AFWL  and  were  test  runs  of  the  FLIR  tracking  a 
receding  aircraft's  jet  exhaust.  The  background  was  clear 
sky.  Other  target  scenarios  with  various  backgrounds  should 
also  be  studied  when  the  data  is  available.  The  four  data 
tapes  supplied  by  AFWL  were  previewed  but  none  were  found  to 
have  anything  other  than  clear  sky  backgrounds. 

REFERENCING  was  used  as  a  tool  to  compute  the  temporal 
cross-correlations  of  a  16-by-16  section  of  sky  between  two 
consecutive  data  frames  (one  thirtieth  of  a  second  apart) . 
REFERENCING  works  by  starting  with  a  16-by-16  reference 
matrix,  Q,  in  the  first  frame.  Assuming  that  the  upper  left 
hand  pixel  coordinates  of  Q  are  (m,n) ,  then  the  cross-cor¬ 
relations  of  Q  with  256  sixteen-by-sixteen  submatrices  of  a 
32-by-32  search  matrix  in  the  next  frame  are  computed. 

Since  the  upper  left-hand  coordinates  of  the  search  matrix 
are  (m-8,  n-8) ,  it  uniformly  surrounds  the  reference  matrix. 
REFERENCING  was  run  so  that  sections  of  clear  sky  background 
FLIR  data  in  one  frame  (i.e.,  specifically  rejecting  target 
intensity  effects)  would  be  searched  for  the  highest  cross 
correlation  to  the  16-by-16  clear  sky  background  reference 
data  of  the  preceding  frame.  The  hypothesis  is  that  the 
same  piece  of  physical  background  imaged  in  Q  has  moved 
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less  than  eight  pixels  (a  good  assumption  based  on  the  tar¬ 
get  scenario)  in  one  thirtieth  of  a  second,  and  that  the 
location  of  this  square  segment  of  background  can  be  found 
in  the  search  matrix  based  on  the  maximum  cross-correlation 
The  formula  for  cross-correlation  coefficient  computations 
is  derived  from  the  standard  statistical  formula  used  to 
determine  correlation  coefficients  and  is  (Ref  13:35) 


Cov[P,Q] 

aPaQ 


E{P,Q}  - 


Vq 


aP°Q 


(C-l 


where  P  is  a  16-by-16  search  submatrix  in  the  frame  follow¬ 
ing  that  of  Q,  and  Cov[P,Q]  is  the  covariance  of  matrices 
P  and  Q  which  is  defined  atypically  in  the  following  way. 
The  expected  value  of  P  and  Q,  E{P,Q} ,  and  their  separate 
means,  nip  and  m^,  are  defined  as  (Ref  12:36): 


16 

16 

E{P,Q) 

256 

£ 

i-1 

£ 

j-1 

Pij 

16 

16 

*P 

£ 

i-1 

£ 

j-1 

Pij 

16 

16 

”Q 

" 

£ 

i-1 

£ 

j-1 

Q|4 


(C-2: 


(c-3: 


(C— 4  ] 


The  definitions  of  E{P,Q},  nip,  mg,  and,  consequently, 
Cov[P,QJ  are  somewhat  unusual  because  they  are  matrices  and 
not  vectors.  crp  and  Oq,  the  standard  deviations  of  P  and  Q 
respectively,  are  computed  as  follows  (Ref  13:36): 
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P 


1 

256 


16  16 
£  £ 
i-1  j-1 


ij 
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The  results  of  the  cross-correlation  runs  are  shown  in 
Table  C-I .  The  columns  of  data  shown  are  the  mission  label 
of  the  data  tape,  the  frames  compared,  the  upper  left-hand 
coordinates  of  the  reference  and  the  search  matrix,  and  the 
maximum  cross-correlation  coefficient.  Table  C-I  shows 
three  separate  runs  as  separated  by  the  dashed  lines.  Each 
run  processed  five  consecutive  frames,  yielding  four  maxi¬ 
mum  cross-correlation  coefficients.  The  starting  upper 
left-hand  coordinates  of  the  first  reference  matrix  in  the 
first  frame  is  a  program  input.  Subsequent  reference 
matrices  in  the  same  run  are  centered  about  the  position  of 
the  maximum  cross-correlation  coefficient.  Thus,  the  upper 
left-hand  coordinates  of  the  reference  matrices  for  the 
second,  third,  and  fourth  computations  of  each  run  may  vary 
depending  on  the  position  of  the  previous  maximum  cross¬ 
correlation  coefficient. 

The  data  shows  a  consistent  positive  trend,  though  it 
is  disappointingly  low,  averaging  .364.  The  clear  sky  back¬ 
ground  has  little  infrared  structure  which  would  remain 
temporally  constant,  and  so  the  low  cross-correlations  are 
not  too  surprising.  More  precise  cross-correlation  coeffi¬ 
cient  calculations  need  to  be  done  on  a  variety  of  fixed 
backgrounds  with  a  motionless  sensor.  The  analysis  per¬ 
formed  establishes  the  temporal  nature  of  the  specific  clear 
sky  background  noise  and  no  more. 
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TABLE  C-I 

Start 

Frames 

Coordinates 

Maximum 

Data  Tape 

Compared 

(Row, 

Column) 

Coefficient 

MIS  10A 

31 

25 

20 

.403 

32 

17 

12 

MIS  10A 

32 

30 

20 

.337 

33 

22 

12 

MIS  10A 

33 

34 

19 

.362 

34 

26 

11 

MIS  10A 

34 

39 

19 

.397 

MIS  10A 
MIS  10A 
MIS  10A 
MIS  10A 


MIS  24A 
MIS  24A 
MIS  24A 
MIS  24A 


- 


1 t 


where  aQ  «  the  rms  value  of  the  dynamics  driving  noise 
At  ■  ti+i  “  fci  “  sample  period 

Using  these  values  is  evaluated  using  the  following 
formula  (Ref  8:171) 


Srd 


/ 


i+1 


£T(ti+i-T)GT(T)QT(T)Gj(T)*J(ti+i-T)dT  (D-4) 


Before  actually  performing  the  integral,  the  matrices  must 
be  multiplied  out.  If  the  following  substitutions  are  made 
in  (0-1) ,  the  multiplications  are  less  confusing: 


A 

B 

C 


e-a(ti+i-ti) 


e"b(ti+l"ti) 


(ti+rti)  e_b(ti+l"ti) 


(D-5) 

(0-6) 

(D-7) 


so  that 
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Using  ♦T(ti+1-ti)  from  CD-8),  from  (D-2)  and  0^  as  shown 
in  CD-3) ,  then 


♦  G  Q 
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| 

1 

1 

1 

mll 

m12 

“is  | 
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1 
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0 

( 

1 

1 
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1 

1 

1 

1 

1 

1 

mn 
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m13 

0 

1 

1 

1 

0 

1 

1 
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m12 

m22 

m23 

1 

1 

1 

t 

1 

1 

m13 

m23 

m33 

(D— 9 ) 


where 


lll 

* 

a2g2 

(D-10) 

l12 

as 

abgiG2 

+  ACG1G3 

(D-ll) 

l13 

as 

ABG1G3 

(D-12) 

l22 

'  m 

B2g2  + 

2BCG2G3  +  C2G2 

(D-13) 

l23 

m‘ 

b2q2G3 

+  BCG2 

(D-14) 

l33 

m 

b2g2 

(D-15) 

If  (D-10)  through  (D-15)  and  CD-5)  through  (D-7)  were  sub¬ 
stituted  into  (D-9) ,  the  resulting  matrix  would  be  cumbersome 
in  size  and  complexity. 

Since  the  integral  of  the  matrix  shown  in  CD-9)  equals 
a  matrix  of  the  integrals  of  each  element,  the  remaining 
demonstration  of  0Td  computation  will  follow  the  integration 
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of  m22  of  (D-9) .  Using  the  expressions  for  B  and  C  from 
CD-6)  and  (D-7) ,  the  expression  for  m22  as  shown  in  (D-13) 
becomes 

m22  “  G|e“2b(ti+l“ti)  +  2G2G3Cti+l”ti)  e“2b  (ti+l"ti) 


+  G2 (ti+1— tA) 2  e"2b<ti+l"ti) 


(D-16) 


The  integration  of  m22  which  would  supply  two  terms  on  the 
main  diagonal  of  is  accomplished  as  follows: 


f'  i+1 

/  *22 lT 
ti 


)dr 


/ti+l 

tG 


2e-2b(ti+1-T)  +  2G2G3(ti+1-x)  e-2b(ti+l”T) 


+  G2(ti+1-x)2  e_2b(ti+l“T)]dx 


(D-17) 


-  G2  J  1  e“2b(ti+l"T)dx+2G2G3yi  1(ti+1-x)e"2b(ti+l“T' 


+  g: 


/^i+i 

(ti+i 


-x)2  e“2b(ti+l“T)dx 


After  computing  the  indicated  integrals  and  regrouping  the 
results,  the  resulting  integration  of  m22  is: 

^i+l 


/i+1 

m22(x)dx 


-  (l-e“2b(ti+l“ti)) [G2/2b  +  2G2G3/(2b)2  +  2G2/(2b)3] 

CD-18) 

+  Cti+1-ti)e"2b(ti+rti)  [G2G3/b+2G2/(2b)2-(ti+1-ti)G2/2b] 


i1* 


dx 
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This  is  the  final  computed  form  of  m22 .  An  equivalent  form 
of  this  equation  can  be  found  in  the  computer  codes  found 
in  Appendix  F.  The  results  of  the  integration  of  the 
remaining  terms  in  (D-9)  can  also  be  found  in  Appendix  F 
and  will  not  be  shown  here. 
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Appendix  E 

Computer  Simulation  Considerations 

To  understand  the  detailed  test  procedure,  a  discussion 
of  the  computer  simulation  techniques  is  appropriate.  The 
software  was  written  in  FORTRAN  IV,  and  the  simulations  were 
accomplished  on  a  CDC  6600.  Due  to  the  size  of  the  software 
and  the  time  .required,  the  software  was  written  (as  origi¬ 
nally  organized  by  Captain  Mercier,  Ref  5:39)  in  three  large 
program  modules  and  was  run  in  two  steps,  to  be  described 
herein.  The  generation  of  Gaussian  noise,  the  approxima¬ 
tions  used  to  produce  FLIR  values,  the  Monte  Carlo  simula¬ 
tion  technique,  and  the  tuning  of  the  filter  are  also  clar¬ 
ified  in  this  discussion. 

For  the  first  step  of  the  computer  simulation,  two 
software  modules,  called  the  main  program  and  the  library 
of  subroutines,  (deck  listings  can  be  found  in  Appendix  F) 
are  run  together.  The  library  of  subroutines  (these  rou¬ 
tines  are  90%  intact  from  their  original  design  by  Mercier, 
Ref  5:85-94)  is  a  group  of  small,  flexible  programs  which 
usually  require  infrequent  modifications.  The  main  program 
contains  the  software  for  setting  up  and  executing  the 
truth  model  and  filter.  Consequently,  this  module  changed 
considerably  and  constantly  during  the  truth  model  evolution 
and  filter  design.  In  its  final  version  the  main  program 
required  about  two  dozen  parameters,  such  as  truth  model  and 
filter  target  sizes  and  intensities,  and  rms  noise  values, 
to  be  read  in  as  data  before  execution.  Besides  requiring 
such  inherent  math  functions  as  sine,  cosine,  absolute 
value,  etc.,  the  main  program  utilizes  a  matrix  inversion 
routine  from  the  IMSL  library  which  is  available  at  AFIT. 

The  library  of  subroutines  also  made  use  of  a  pseudo-random 
number  generator. 

The  third  module,  the  PLOT  routines,  (the  PLOT  deck 
listing  is  found  in  Appendix  G)  is  run  subsequent  to  the 
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execution  of  the  main  program  and  the  library  of  subrou¬ 
tines,  which  generate  and  store  the  data  used  by  the  PLOT 
routines.  These  programs,  which  were  written  by  Captain 
Mercier  in  1978,  required  only  small  modifications  for  their 
usage  in  this  study  to  generate  the  plots  seen  in  this 
Appendix,  in  Chapters  III,  V,  and  VI,  and  in  Appendices  I, 

J,  and  L. 

The  Gaussian  white  noise  terms  used  in  the  simulation 
are  generated  in  a  subroutine  using  a  pseudo-random  number 
generator  with  a  uniform  probability  density  function 
between  0  and  1.  To  produce  a  discrete  realization,  v(k) 
of  Gaussian  distributed,  white  noise  of  a  given  variance, 
ow,  twelve  independent  values  are  used  from  the  random  num¬ 
ber  generator  in  the  following  equation  (Ref  5:41): 

12 

v(k)  *  a  {  Z  u[12 (k-1)  +  i]  -6}  (E-l) 

W  i=l 

for  k  —  1,2,3,... 

where  u(k)  =  random  number  from  the  RAND  function. 

This  equation  will  produce  a  white  sequence  of  numbers  that 
is  very  nearly  Gaussian  with  zero  mean  and  desired  variance, 
as  can  be  argued  based  on  the  Central  Limit  Theorem  (Ref  8: 
109-110) . 


Reasonable  approximation  methods  for  evaluating  the 
outputs  of  the  FLIR  are  needed  for  the  computer  simulation. 
The  precise  mathematical  equation  would  entail  evaluating 


27 


/ 


7* 


as  in  equation  (7)  for  the  truth  model.  Similar  expressions 
are  evaluated  in  the  filter's  estimate  of  the  measurements 
and  in  the  computation  of  the  partials  found  in  equation 
(15) .  Equation  (E-2)  above  represents  the  average  value  of 
the  Gaussian  intensity  fraction  for  the  pixel.  For  online 
usage  a  reasonable  approximation  of  equation  (E-2)  is  simply 
to  evaluate  the  intensity  function 


(E-4) 


at  the  middle  of  the  pixel,  and  the  filter  uses  this  single 
point  approximation  to  the  integral.  This  approximation  is 
poor  when  the  eigenvalues  of  P  are  small  compared  to  the 
size  of  a  pixel.  A  more  accurate  approximation  can  be  com¬ 
puted  by  dividing  the  pixel  into  smaller  areas,  computing 
the  midpoint  value  for  each  of  these  small  areas,  and  aver¬ 
aging  the  results.  Because  the  truth  model  is  to  depict 
the  real  world,  an  accurate  pixel  evaluation  was  performed 
by  using  the  midpoints  of  a  4-by-4  array  of  equal  subdivi¬ 
sions  of  each  pixel. 

As  argued  extensively  by  Captain  Mercier  on  page  41  of 
his  thesis,  a  Monte  Carlo  simulation  is  needed  to  accomplish 
a  performance  evaluation.  The  inherent  goal  of  Monte  Carlo 
simulation  technique  is  that  sufficient  numbers  of  simula¬ 
tions  be  run  so  that  convergence  is  reached  in  the  standard 
deviation  (and  other  estimated  statistics)  of  filter  errors, 
as  shown  in  Figure  E-l.  A  simple  covariance  analysis  of 
filter  performance  is  not  possible  for  the  adaptive  extended 
Kalman  filter  investigated  in  this  report  due  to  the  truth 
model  nonlinearities  involved  (Ref  8:329).  Consequently, 
a  sample-by-sample  simulation  using  white  Gaussian  noise 
realizations  generated  by  (E-l)  is  conducted.  The  adaptive 
extended  Kalman  filter  produces  state  estimates  from  the 
available  noisy  measurements.  If,  after  n  number  of 
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simulations,  the  standard  deviation  of  the  filter  errors 
and  other  estimated  statistics  show  convergence,  a  meaning¬ 
ful  performance  evaluation  can  be  achieved,  and  the  Monte 
Carlo  simulation  can  be  terminated.  Captain  Mercier  per¬ 
formed  such  runs  and  found  that  20  (n)  simulations  were 
sufficient  to  achieve  such  convergence  (Ref  5:42)  as  demon¬ 
strated  in  Figure  E-l. 

Tuning,  which  is  a  process  performed  off  line  by  the 
designer  is  aimed  at  selecting  the  correlations  times  and 
noise  strengths  to  allow  effective  filter  performance. 
Effective  robust  tuning  will  allow  the  filter  to  perform 
acceptable  over  a  range  of  possible  design  conditions.  Cap¬ 
tain  Mercier  performed  this  tuning  by  adjusting  Qpd  and  R^ 
to  minimize  state  errors  in  an  appropriate  manner.  Depend¬ 
ing  on  the  problem  the  objective  may  be  to  minimize  the 
error  variance  of  all  states  simultaneously  or  to  minimize 
only  a  given  number  of  states'  error  variances  (Ref  5:44). 

(  The  proper  tuning  of  the  filter  is  indicated  by  agreement 

between  the  filter's  estimate  of  the  variance  of  the  error 
with  the  actual  error  variance  for  each  state.  Case  2  in 
Appendix  I,  which  is  the  baseline  case  for  the  robustness 
runs  in  Chapter  III,  demonstrates  a  properly  tuned  filter. 

By  the  above  definitions,  the  adaptive  estimation  of 
formulated  in  Chapter  IV  and  discussed  in  Chapters  V  and  VI 
can  be  viewed  as  a  robust  self-tuning  technique. 
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This  appendix  contains  the  FORTRAN  code  for  the  algo¬ 
rithms  presented  earlier  in  the  text.  Presented  first  is 
the  eight  state  filter  program  in  its  final  form.  Next, 
the  six  state  filter  program  in  its  final  form  is  listed. 
The  third  listing  is  of  the  72-state  temporal  background 
noise,  truth  model  with  Mercier's  four  state  filter.  Fin¬ 
ally,  the  subroutines  used  by  the  eight  state  program  are 
listed.  The  eight  state  filter  subroutines  cannot  be  used 
with  the  six  or  four  state  filter  as  minor  modifications  to 
several  subroutines  are  required  for  use  with  the  six  or 
four  state  filter. 
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4 

# 

4 

0 


,8)  ,  SAVE(l.A),QPDi(8,8), 

,8),XS(8)  ,  H  (6L  ,  8 )  ,  PFP  (8,8)  , 
,RIH(6U,1),  OXDXT (8,8), PPFP (8,  6)  » 
*>4>  , RN(84,6A)  ,80(8,2)  ,PHIFT(8,8) 
!2) ,PFP0LDCe,8) tXFP0(8) 


RE  AC  ANO  ECHO  DATA 


REAO  *»  SIGS1 


PRINT  «,  "RMS  OYHAHICS  FOR  TRUTH  HOOEL, 
REAO  *,  SIGHA3 

SIGS1 

S 

", SIGS1 

PRINT  *,  "RMS  TRUTH  MOOEL  BACKGROUND  NOI 
REAO  *,  SIGFLR 

SE, 

SIGMAB 

s 

", SIGMAB 

PRINT  ♦,  “RMS  TRUTH  HOPEL  FLIR  N0ISEy 
REAO  ♦,  IMAX 

SIGFLR 

s 

“,SIGFLR 

PRINT  ♦,  "TRUTH  MOOEL  MAX  INTENSITY, 

READ  *,  SIGAT 

IMAX 

c 

", IMAX 

PRINT  *,  “RMS  ATMOSPHERICS  FOR  TRUTH  MOOEL, 
REAO  *,  NRUN 

SIGAT 

s 

", SIGAT 

PRINT  *,  "NUMBER  OF  MONTE  CARLO  RUNS, 
REAO  ♦,  T FINAL 

NRUN 

s 

",NRUN 

PRINT  *,  "FINAL  TIME, 

REAO  «, SIGHS 

TFINAL 

s 

", TFINAL 

PRINT  *,  "INITIAL  RMS  TRUTH  MODEw  SIGHA 
REAO  *,  ASPRO 

PERVEL 

,  SIGHS 

s 

", SIGHS 

PRINT  ♦,  "TARGET  AS®ECT  RATIO, 

REAO  *,X0 

ASPRO 

s 

", ASPRO 

PRINT  *,  "INITIAL  X  POSITION, 

REAO  * ,Y0 

xo 

s 

",X0 

PRINT  *,  “INITIAL  Y  POSITION, 

REAO  *,Z0 

YO 

s 

",  YO 

PRINT  *,  "INITIAL  ?  POSITION, 

REAO  *, XOOTO 

ZO 

s 

",Z0 

PRINT  ♦,  "INITIAL  VELOCITY,  X  DIRECTION 
REAO  *, YOOTO 

• 

XOOTO 

8 

", XOOTO 

PRINT  *,  "INITIAL  VELOCITY,  Y  DIRECTION 
REAO  ♦,ZOOTO 

• 

YOOTO 

8 

", YOOTO 

PRINT  ♦,  "INITIAL  VELOCITY,  Z  DIRECTION 
REAO  ♦, ISPTL 

• 

ZDOTO 

8 

",ZOOT  0 

PRINT  *,  "SPATIAL  NOISE*  1-YES,  8-NO 

ISPTL 

8 

", ISPTL 

( 


IF(ISPTL.NE.I)  GO  TO  2 
REAO  *,C(l),C(2),C(3),r(4),C(5) 

PRINT  ♦,  "SPATIAL  NOISE  CORRELATION  COEFFICIENTS*  " 

PRINT  +  ,C(1),C(2)  yC(3)  ,C(4),C(5) 

CONTI  NUE 
PRINT  ♦,"  " 

REAO  *,  SIGHFO 

PRINT  *,  "INITIAL  FILTFR  SIGMA  VELOCITY,  SIGHFO  « 
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RE AO  *»  ARO 

PRINT  *,  -INITIAL  FILTER  ASPECT  RATIO,  '  ARO  *  -,ARO 

REAO  *,  SIGF2 

PRINT  *,  “RHS  FILTER  ATMOSPHERIC  NOISE,  SICF2  *  -,SIGF2 

REAO  ♦,  RF 

PRINT  ♦,  -FILTER  MEASUREMENT  NOISE  RMS,  RF  *  -,RF 

REAO  ♦,  SIGF10 

PRINT  *,  -INITIAL  RMS  OYNAMICS  FOR  FILTER,  SIGF1&  *  -,SIGF10 
REAO  ♦,  FIMAXO 

PRINT  *,  "FILTER  MAX  INTENSITY,  FIMAXO  *  -, FIMAXO 

REAO  *,IPRINT 

PRINT  *,  -PRINT  COMMAND,  i  YES,  0  NO  IPRINT  *  -,IPRINT 

PARAMETER  VALUES 

YF0V*8. 

NPIX*8 
ATAU1*14.14 
AT  AU2-659  .5 
FTAU2*1./ATAU1 
OT  =  1./30. 

IREF*1 

FI MAX*ABS (FIMAXO) 

SI  GVF * ABS  (SI  GMF 1) 

SIGF1*A8S (SIGF13) 

QFOMAX  (1)  *2. 

OFOMA  X (2) *QFOMAX(l) 

0F0MAX<3)*14. 

QFOMAX (4)  *QFDNAX(  3) 

OFOMA  X  (51  *20. 

QFOMAX  (6)  *QFOMAX  (5) 

QFOMAX (7)  **5 
QFOMAX (6)  *QF0MAX(7) 

INITIALIZE  TRUTH  MODEL  VARIA9.ES 

CALL  RANSET(75632) 

ONE  *  1 
NPS  *  8 
NMS  *  NPIX**2 
NFS  *8 

NFS2*NFS*NFS 

NFSM2*NFS-2 

NIS  *  3 

RI*1./RF 

SN*SIGMA8/IMAX 

IF(SN.LE.O.)  SN*.  031 

SN*1,  /SN 

RANGE  0*SQRT  CX3**2+Y0**?*Z0**2) 

AGAIN  *  .351036934  •  STGAT 

IFILE-6 

OELT  *  -l.*OT 

DO  5  1*1,8 

BO (I, 1)*0 • 

BO (I, 2)*0 • 

00  5  J*l,  8 
00  (I,  J)  *  0. 


SQQOdf J)  a  0. 

PHI  (Z  f  J)  *  0. 

CONTINUE 

FACT*(AGAIN**2)*(ATA'J1**2)*(ATAU2**'4) 

FACT1  *  ATAU1-ATAU2 

FACT2  a  ATAU1+A TAU2 

FACTS  a  2.*ATAU2 

XF0V*8, 

61  a  FACT/ (FACT 1**4) 

62  «  FACT/(FA5T 1**3) 

G3  a  FACT/CFACT 1**2) 

R1  *  1.-  EXP(2.*ATAUl'OELT) 

R2  *  1.-  EXP(FACT2*0ELT) 

R3  a  l.-  EXP(2 •* ATA U2*  CELT) 

R4  *  OT* EXP (CELT* FAC T2) 

R5  *  CT*EXP(2.*ATAU2*0FLT) 

FILL  OUT  TRUTH  MOOEL  PHI  MATRIX. 

SEE  MERCIER'S  THESIS  FOR  DERIVATION 


PHI  (1,1): 
PHI <2  »  2) 
PHI  (3  ,  3) 
PHI  (4,4) 
PH  1(4  ,5) 
PHI  (E  ,5) 
PHI  (6, 6) 
PHI  (7  ,7) 
PHI  (7  ,8) 
PHI  (8  ,8) 


1. 

‘  PHI  (1,1) 

EXP ( AT  AUl*OFLT) 
EX® ( AT  AU2*QFLT) 
OELT*PHI(4,t ) 
PHI (4,4) 

PHI (3, 3) 

PHI (4, 4) 

PHI (4, 5) 

PHI (5,5) 


WRITE  (6,11) 

L  FORMAT (///2X,"THE  TRUTH  MODEL  STATE  TRANSITION  MATRIX  IS*"/) 
CALL  MCUT (PHI, NPS,NPS) 

FILL  OUT  OISCRETE  INPUT  MATRIX 

BO  (1,1) «  OT 
80(2,  2)=  OT 
WRITE (6,15) 

FORMAT (///2X, "THE  TRUTH  MOOEL  INPUT  MATRIX  IS*"/) 

CALL  M0UT(B0,NPS,2) 

FILL  THE  QO  MATRIX  WITH  VALUES  USING  EXACT  INTEGRATION 
SEE  MERCIER'S  THESIS  FOR  DERIVATION 

00(1,1)*  SIGS1 
QO  (2,  2)  a  Q0(1,1) 

00(3,3)  *  (G1*R1)/(2.*ATAU1) 

00(3,4)  «  R2*(G2/FACT2**2-G1/FACT2)-P4*G2/FACT2 
00(3,5)  *  G2*R2/F ACT 2 
00(4,3)  a  00(3,4) 

00(4,4)  *  R3*(Gi/FACT3-2,*G2/FACT3**2*2.*G3/FACT3**3)- 
*  R5*(G2/ATAU24G3*0T/FACT3-2.*G3/?ACT3**2) 

00(4,5)  a  R3*(G3/FACT3**2-G2/FACf  3)  «R5* 63/ FACTS 
00(5,3)  «  00(3,5) 

00(5,4)  >  00(4,5) 

00(5,5)  «  R3*G3 /FACT  3 
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I  00  20  1*3,5  . 

00  2C  J*3 ,5 
00  (1  +  3,  J*3)  *  00(1,  J) 

0CI-2, J-2)  *  00(1, J) 

28  CONTINUE 
WRITE (6, 30) 

'  30  FORMA T(///2X, -THE  TRUTH  MODEL  QD  MATRIX  IS*"/) 

CALL  MOUT  (  Q0,NPS,NPS) 

TAKING  CHOLESKY  SQUARE  ROOT  0r  00 

SQ00( 1,1)  *  SORT (Q0( 1,1) ) 

SQOO(  2,2)  *  SQQ0(  1,  1 ) 

CALL  CH0LESK(Q,  V0RK,NI?) 

00  33  1*1 ,NIS 
00  33  J*1  ,NIS 
SQQ0( I+2»  J+2)  *  WORK  (J,I) 

SQQ0(I+5, J+5)  *  WORK ( J,I) 

33  CONTINUE 
WRITE (6,35) 

35  FORMAT (///2X, -THE  CHOLFSKY  SQUARE  ROOT  OF  00  IS*"/) 
CALL  MOUT (SQQO,  NPS,  N°S ) 

IF  (ISPTL.NE.l)  GO  TO  *1 

SET  UP  SPATIAL  NOISE  CORRELATION  COEFFICIENT"  MATRIX 

N*64 
M*8 

00  36  1*1, N 
R(I,I )*1« 

IF  (I.GE.64)  GO  TO  36 
R(I,I+1)*C(1) 

IF  (I.GE.63)  GO  TO  36 
R(I,I +2)  =C  (3) 

IF  (I.GE.57)  GO  TO  36 
R( 1,1 *61*0(4) 

R(  I, I  +7)  *C(2) 

R(I,I+6)*C(1) 

R( I, I +9) *C (2) 

R(  I, I +  1*1)  *C(4) 

IF  (I.GE.49)  SO  TO  36 
R( 1,1+14) *C(5) 

R( 1,1+15) *C(4) 

R(  1,1  +16)  *C(3) 

R(  1,1+17)  «C<4> 

R(I,I+lfl) *C(5) 

36  CONTINUE 

00  37  1*1  ,M 
RC8+I -7,e+I)*8.0 
RC8+I-7,8+I-l>*0.0 
R( 6+1 -6,8+1) *0.0 
IF  (I.GE.8)  GO  TO  37 
R( 8+1 ,8+1 +1) *0*8 
RC 8*1, 8*1 +2) *0. 3 
R(6+X-l,8+X+l>*0.9 
R(G+X-7,8+I+7) *0*0 
R<8+I-7,8+I+8)*0.0 
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FTN  >.7*476 


RC8*I -6,8*I*8)  0.0 
IF  (1.GE.7)  GO  TO  37 
(  RC  6*1 f  8*1 *9) *0.3 

R< 8*1 ,8*1 *10) *0.0 
R<  8*1  -1*8*1  *9)  0 

ZF  (Z.GE.6)  GO  TO  37 
RC8*I  ,8*1  *17) *0 .C 
RC  8*1  ,8*1 *18) *0.6 
RC8*I-l,8*I*17)«C.O 

37  .  CONTINUE 

DO  38  1*1, N 
L«I*1 

00  38  J*L  »N 
RCJ,I)*RCI,J) 

38  CONTZNUE 

00  39  Z*1  »N 
DO  39  J*i»N 
RNCZ,  J)*SIGMA3*RCI,  J) 

39  CONTZNUE 
MRZTE (6,44) 

44  FORMAT C///2X,-64  X  64  SPATIAL  CORRELATION  MATRIX!-/) 
CALL  M WRITE  CRN, 84, 64, (4) 

C 

C  COMPUTE  THE  CH  OLE  SKY  SOUARE  ROOT  OF  RN 

C 

CALL  CHOLESK  CRN,  P.,64 ) 

WRITE (6*48) 

(  48  FORMAT (///2X, "THE  CHOLESKY  SOUAR*  ROOT  OF  RNt-/) 

CALL  M WRITE  CR,64,64,6M 
GO  TO  42 

41  00  43  1*1*64 

43  R(I»I)*1. 

42  CONTZNUE 
C 

C  SET  UP  FILTER  MATRICES. 

C 

DO  51  I*1*NFS 
00  51  J*i  *NFS 
PHIFd*  J)  *0. 

51  QFOCI ,  J)*0. 

C 

C  FILL  OUT  FILTER  PHI  MATRIX 

C 

PHIFC 1*1) *1. 

PMIF< 1*3) *0T 
•  PHIFC  1,5) *0T**2/2. 

PHIFC 2,2) *PHZF (1,1) 

PHIFC 2,4) *0T 
PHIFC 2*6)  «PHIFC  1,5) 

PHIFC  3,3) *PHIFC1, 1) 

PHIFC 9,5) *0T 
PHIFC4.4) *PHIF (1,1) 

PHIFC4,6>  *0T 
PHIFC5,5) *1* 

™  PHIFC6,6)*1. 

PHIFC 7 ,7) *EXPC0ELT/CTAU2) 

PHIFC  8,8)  *PHIF  C7,7) 
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00  56  I*i,NFS 
00  56  J*i  »NFS 
PHIFT(I,J)*PHIF(J,I) 

FILL  OUT  FILTER  DISCRETE  OFD  NATRIX 
FOR  START  OF  ACQUISITION  PHASE 


0FD(i  ,i)«0T**5*SIGFi/2P. 

QF0(l,3)»0T**4*SIGFi/8. 

0F0(i,5)*DT*»3*$IGFl/6. 

0F0C2,2)*0F0(i,l> 

QF0(2  ,4)*QFD (i,  3) 

OF 0(2  ,6)  =  QF0  (It  5) 

0F0(3,i>*0FD(l,3) 

OF 0(3 , 3) *0T**3*  SIGFi /3 . 

0F0(3,5)=0T**2*SIGFi/2. 

QFO  (4  ,2)* QFO (3,  1) 

QF0(4 ,4) =QFD (3, 3) 

0F0(4,6>*QFD(3,5) 

0F0(5 ,1)  =  QF0(1,5) 

OF 0(5  »3)=QFD  C3»5) 

0F0(5,5)  =  SIGFi*0T 
QFO  (6  »  2) * QFO (1,5) 

OFO  (6 , 4)  a QFO (3*5) 

0F0(6,6)sQFD(5,5) 

QFO  (7,7)  *  (SIGF2**2>'(1*-EXP(2.*0ELT/FT#U2)) 

QFO  (8,8) *QFD  (7,7) 

WRITE (6,40) 

)  FORMAT (///2X, "THE  FILTFR  STATE  TRANSITION  MATRIX  IS*"/) 

CALL  MOUT (PHIF, NFS,  N  PS) 

WRITE (6,45) 

45  FORMAT (///2X, "THE  INITIAL  FILTER  QO  MATRIX  IS*"/) 

CALL  MOUT (  QFO, NFS, NFS) 

PRINT  *,"  " 

PRINT  " 


PRINT  ♦,-#**»*****»  9EGIN  THE  MONTE  CARLO  SIMULATION 

00  99  L*1,NRUN 
TIME  *  0. 

XCENTRaO. 

YCENTRaO# 

FI MIN  «0. 

IMEASaQ 

FI MAX*A6S (FIMAX  8) 

ARsARC 

SI GVF  * ABS (SIGMF  0) 

SI GF1*ABS (SIGF13) 

VARYQaABS (SIGFi fl) 

TRXXT *200  • 

KANIN  O«0 


RESET  INITIAL  CONDITIONS  FOR  NEW  RUN 

DO  46  I*i,  NPS 
XS(I)  *  0. 

46  CONTINUE 
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00  47  1*1 .NFS  -  • 

XFP(I )  *0. 

XFM(X)*0. 

DC  47  J*i,NFS 
QF0(I ,  J)*0. 

PFP(X.J)  *  0. 

47  CONTINUE 

PROVIDE  FILTER  WITH  INITIAL  VELOCITY  AND  POSITION 

RHOR*  CX3**2+73**2) 

RANGF*(RHOR*Y3**2) 

XFP(3)*(Z0*XOOT0-X5*7OCT0)/(RHOR*.000  0  2) 

RH OR* SORT  (RHOR) 

XFP(4  )*(RHOR*YOOTO-Y  0*(  (XO*XOOTOfZO*ZDOTO)  /RHOR) )  /  (RANGE*.  uOO  02) 
IFCIFRINT.NE.l)  GO  TO  1020 

PRINT  RHOR*", RHOR,"  RANGE**2»", RANGE,"  XFP(3)*",XFP(3) , 

*"  XFP(4)  *",XFP  (4) 

1020  CONTINUE 

FILL  IN  P*  AT  TINS  P 


PFP(1,1)  =  25. 

PFP(2,2)*PFP(1, 1) 

PFP(3 ,3)*2000« 

PFP<4 ,4)*PFP(3,  3) 

PFPC?,5)*100. 

PFP(6 ,6)*100. 

PFPC7  ,7)*.2 
PFP(8,S)*«2 

FILL  IN  QFO  AT  TINE  0 

QFOd  , 1)  *0T**5*SIGF1/2P. 

QFOC1 , 3)*0T**4*SI  GF1/8  • 

OFOd  ,5) *OT**3*SIGFl/6. 

OFO(2,2)*OFO  (1, 1) 

QFO (2 ,4) * QFO (1, 3) 

QFO (2 , 6)* QFO (1, 5) 

OF 0(3 ,1)*0F0(1, 3) 

0F0(3, 3) *0T**3*SIGFi/3. 

QFO (3 ,5)*0T**2*SIGFl/2. 

QFO  (4 ,2)* QFO  (3,1) 

QF0(4,4)*QF0(3,  3) 

OF 0(4,6)* QFO (3, 5) 

QFO  (5 ,1)*0FD(1,  5) 

OF  0(5 , 3)  *  QFO  (3, 5) 

QFD(5  ,5)*SIGF1*0T 
OF 0(6 ,2)*  QFO  (1,5) 

QFO  (6 ,4)*QF0(3,  5) 

QFO (6, 6)* QFO (5,  5) 

QFO  (7,7)  *  (SIGF2**2)M1.-£XP(2.*0ELT/FTAU2)) 
QFO  ( 8,8)*QF0  (7,7) 

C 
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TIME  LOOP  STARTS  HERE 

$0  TIME  *  TIME  ♦  OT 

IF(TIME.GT.TFINAL)  SO  TO  99 

PERFORM  TR'JTH  MOTEL  SIMULATION 


VTIHE*TIME-0T/2. 

XV EH* -500 .*TIME*XQ 
YVEH*-300.*TI1E»Y0 
IV EH* 20 
XOOT* -500 • 

TOOT* -300. 

70  OT*  0  • 

RHOR* (XVEH**2+ZVEH‘*  2) 

RANGE*  (RH0R*YVEH**2) 

UT  (1,  1)*(Z.VEH*XOOT-XVEH*ZOOT)  /(RHOR*.  00002) 

RHOR* SORT  (RHOR) 

UT  (2, l)*(RHOR*YOOT-Y VEH* ( (XVEH*X3 OT+Z  VEH*ZDOT) /RHOR) ) / (RANGE 
-*.000  02) 

RANGE*SQRT  (RANGE) 

VMAX*SQRT ( XOOT** 2 ♦Y9 OT* *  2+Z00T**2) /(RANGE* .00002) 

IF (IPRINT .NE.l)  GO  TO  1025 
PRINT  *,“  “ 

PRINT  *,“  TIME  *  M,TIHE 

PRINT  *,"XDOT*",XDOT,"  VDOT*",Y30T, -  XVEH*",XVEH,~  YVEH*", YVEH, 
*"  RHOR*", RHOR,"  RANGF*"»RANGE 
PRINT  *,“UT(1,1)«",UT(1,1),"  UT( 2,1) *",UT (2,1) , **  VMAX*",VMAX 
10  25  CONTINUE 

XC  ENT  R*XC  ENTRHIT  ( 1 , 1  )*  BD  ( 1, 1 ) 

YCENT  R*YCENTP>UT(2»  1  )*BO  (2f  2) 

CALL  NOISE  (N»S,W) 

CALL  HMPY(TEHPf SQQOf MfNPSfNPSfONE) 

CALL  HMPY (TEMPI, PHI,  XS,NPS.NPS, ONE) 

CALL  M  AOO  (XS  »  TEMP  fTEf’PlfNPS,  ONE,  ONE) 

CALL  MMPY (TEMPI, PD, UT,NPS»2f ONE) 

CALL  MADO(XS, XS, TEMPI, NPS, ONE, ONE) 

CALL  SHIFTA(XFP,XFPO,NFS, NFS) 

FILTER  STATE  PROPAGATION. 

121  CONTINUE 

CALL  KHPY(XFM,PHIF,XFt>,NFS,NFS,ONE) 

IF (IPRINT .NE.l)  GO  TO  1001 
IF (TIME. LT. 3.0)  GO  TO  1001 
PRINT  *,  "  “ 

PRINT  *,"XS" 

CALL  MOUT (XS,NPS, ONE) 

PRINT  *,  *  * 

PRINT  *,"XFM" 

CALL  MOUT (XFM,NFS, ONE) 

PRINT  *,"  " 

PRINT  *,"  QFO  EST" 

CALL  MCUT (QFO, NFS, NFS) 

1001  CONTINUE 


FORM  CENTROID  POSITION  ANO  FILL  TRUTH  ARRAY 
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XPEAK  *  XS(1)  ♦  XS(3)  ♦  XSC4)  -<FH(1) 

YPEAK  *  XSC2)  ♦  XS(S>  ♦  XS(7)  -<FM(2) 

IF (IFRINT  «N£.l)  GO  TO  1630 

PRINT  *, -XPEAK*-, XPEAK, “  YPEAK*-, YPEAK,"  FIMAX*“,FIMAX 
10  30  CONTINUE 

IF (AES (XPEAK) «GT«  3»*ASPR0#SIGMS J  GO  TO  101 
IF(ABS(YPEAK).GT.3.*ASeR0*SIGMS)  GO  TO  1C2 
CALL  PEAS (XPEAK, YPEAK, A ,Z,R> 

IF (IPRINT  *NE*1)  GO  TO  1002 
IFUIHE.LT.3.8)  GO  TO  1002 
PRINT  *,-  - 
PRINT  *,-Z" 

CALL  MOUT  (Z,N°IX,NPIX) 

1002  CONTINUE 
C 

C  SEARCH  FOR  FIMAX  AND  FIHIN 

C 

IF (FI MAXO •GT«0« )  GO  TO  59 

FIMXOsFIMAX 

FI MAX *0* 

FIHINsC. 

00  55  1*1,  NHS 

IF (Z( I) •GT.FINAX)  FI  t*AX  *  Z(I) 

55  IF  (Z( I ) •LT«FIMIN)  FI  MIN  *  Z(I> 

FIHAX*FIMAX*.14*SIGVF*»2-1.52*SI3VF+U.3«; 

FI MAX*  .8*FIHX04  .2*FIMAX 
59  CONTINUE 
C 

C  FILTER  COVARIANCE  PROPAGATION 

C 

CALL  SHIFTA (PFP, PFPOLD, NFS2,  NFS2) 

CALL  HMPY (EXTRA , P  HI F, PFP , NFS , NFS,  NFS) 

CALL  MMPY (PPFP, EXTRA ,PHIFT, NFS, Nr S, NFS) 

C  PPFP*  PHI  *  P(TI -1) ♦  ♦  PHIT 

CALL  HAOD  (PFM,  PPFP,  0*0,  NFS,  NFS,  OV  E) 

C  PFH  *  P(TI-)*PHI  *  F (TI*i)  ♦  *  PHIT  ♦  G  *  Q  *  GT 

IF (IPRINT «NE«1)  GO  TO  1008 
IF (TI ME.LT.3.8)  GO  TO  1008 
PRINT  ♦,-  - 
PRINT  *,-  PFH- 
CALL  HOUT (PFN,NFS,NFS) 

1008  CONTINUE 
C 

C  PERFORM  MEASUREMENT  UPOATE  FO*  THE  FILTER 

C  INVERSE  COVARIANCE  FORM 

C  FORM  FILTER  CENTROID  POSITION  ANO  FILL  OUT  NON  LINEAR 

C  SMALL  H.  CALCULATE  PARTIAL  SMI LL  H  PARTIAL  X 

C 

XPEAK  «  XFM(7) 

YPEAK  *  X FM ( 8 ) 

IFC(XFM(3)»*2*XFM(4)**2)  .EQ.O.)  <FM(3>*.001 
CALL  MEASF(XPEAK, YPE  AK,ONE,HF,H) 

IF (IPRINT .NE.1)  GO  TO  1003 
IFCTIME.LT.3.S)  GO  TO  1003 
PRINT  ♦,-  - 
PRINT  *,-HF- 


ooo 


CALL  HCUT(HF,MPlX,Nf»IX> 

IOC 3  CONTINUE 

00  60  I*1,NMS 
00  6C  J*1,NFS 
HT (J|I)  * 

60  CONTINUE 

CALL  HMPYCPFP,HT,H,NFS,NMS,NFS) 

C  PFP*  HT  *  H  (R-l  IS  SCALAR  A*0  MULTIPLIED  LATER) 

IOGT  *  0 

CA  LL  L INV 2F(  PFM , NFS , NFS , EXTRA , 103  T , WK ARE  A , IER) 

C  EXTRA*  P(TI-) -1 

00  65  1*1  (NFS 
00  65  J=1 »NFS 

PFP  (I  , J)*PFP(I,  J)*RI  ♦  EXTRA  (It  J) 

65  CONTINUE 

C  PFP*  P (TI*> -1=  HT  *  R-l  *  H  ♦  P(TI-) -1 

IOGT  *  0 

CALL  L INV 2F( PFP t NFS t  NPS»  EXTRA , IDS T,  WK ARE A, IER) 

DO  7C  1*1 1  NFS 
00  70  J*1,NFS 
PFP(I,J)  *  EXTRA(I,J) 

70  CONTINUE 

00  13 C  1= 1»NF3 

1F(PFP(I,I).GT.3.)  GO  TO  130 

PRINT  *,"PFP ("»I, I,")*", PFP  (1,1),"  RUN=~,L,~  TIME*", 
PFP(I,I)  =  PFPOLO(I,I) 

130  CONTINUE 

IF  (IPRINT.NE.l)  GO  TO  1004 
IF(TIME.LT.3.3)  GO  TO  1004 
PRINT  " 

PRINT  *,"PFP" 

CALL  HOUT (PFP, NFS, NFS) 

10  04  CONTINUE 
C  PFP*P(TI*) 

CALL  MAOD (RIH ,Z,HF,NHS,ONE,NFS) 

C  RIH*RESIOUALS*  Z-SMALL  H 

IF  (SIGHFO  .GT ••]•)  GO  TO  62 

COMPUTE  NEW  ESTIMATE  OF  SIGVF  ANO  AR 

00  50  1*1, NMS 

AR*  AR4.m*CRIH<I)*PL?P(I,i)) 

SIGVF *SIGVF+. 3C1*  (RIH( I)  #PL2P(I»2) ) 

50  CONTINUE 

IF(AR.LE.l.)  AR*1. 

IF  (SIGVF. LE.O.)  SIGVFsSIGMFO 
62  CONTINUE 

CALL  MHPY (EXTRA ,HT, RIH, NFS, NMS, 04 E) 

IF  (IPRINT.NE.l)  GO  TO  1007 
IF (TIME.LT.3.S)  GO  TO  1007 
PRINT  " 

PRINT  ♦,"HT  X  (7  -  SMALL  H)" 

CALL  MGUT(  EXTRA,  NFS,  ONE) 

1007  CONTINUE 

C  EXTRA*  HT  *  (7-SMALL  H) 

CALL  NHPY(OX,PFP,  EXTRA, NFS, NFS, OWE) 

C  OX  *  P (TX+)  *  HT  *  (Z  -  SMALL  H> 


TIME 


( 


00  7?  1*1  .NFS 
OX  (I)*DX  CI)*RI 
75  CONTINUE 

C  OX  *OSLTA  X  *  P(TI*)  ♦  HT  ♦  R-l  *  CZ  -  SHALL  H) 

DO  78  1*1, NFS 
00  78  J*1 .NFS 
0X0 XT (I,J)*0X  (I)*DX  (J) 

78  CONTINUE 

IF (IPRINT  .NE.1)  GO  TO  1305 
IF (TI HE.LT .3.8)  GO  TO  1005 
PRINT  “ 

PRINT  *,**  OXOXT" 

CALL  HOUT (OXOXT, NFS, NF?) 

10  05  CONTINUE 

C  OXOXT*  (X (TI ♦)  -  X (TI-))  *  CX(TI*)  -  X(TI«))T 

TRXXT  0*TRXXT 
TRXXT  *C, 

00  79  1=1,  NFS 

79  TRXXT*TRXXT+DXDXT  (1,1) 

IF (TRXXT .LE.13.*T9XXT0)  GO  TO  112 
IF (HA NINO . EQ. 1)  GO  TO  112 
ECC=SQRT(1.-(1./AR**2)) 

SIG(i)=SQRT(SIGPVF**2/(l.-(ECC*C3TH)**2) ) 
SIG(2)=SQRT(SI&PVF**2/(1.-(ECC*SUH>**2)  > 
IFCIPRINT.NE.l)  GO  TO  1040 

PRINT  *,*•  SIG(l)  *n,SIG(l)  ,**  SIG(2)*M,SIG(2)  ,**ECC***,ECC 

1040  CONTINUE 

00  116  1*1,2 
OHTZ*EXTRA(I) 

AL 100*  (ALOGIO  (A  3S (9HT7) )  ♦•4651*SIG(I)  -3.9)/.82 
0(1)* EXP ( 2*303* AL 100) 

IFCIPRINT.NE.l)  GO  TO  1045 
PRINT  ♦, -0(“, I,")***,  0(1),-  OHTZ*  -,0HTZ 
1045  CONTINUE 

XFP  (I  ♦4)=(2.*0(I)/OT**2)*OHTZ/ARS(OHTZ)*XFPO(I+4) 

116  CONTINUE 

00  114  I* l.NFS 
IF (I. EQ.5 )  GO  TO  114 
IF (I. EC. 6)  GO  TO  114 
XFP(I)*XFPO(I) 

114  CONTINUE 

CALL  SHIFTA(PFP0L0,°FP,NFS2,NFS2» 

00  140  1*3,4 

VFACT 0R=10./SQRT(PFP  (1,1)) 

AF ACTOR*300. /SORT (PFP (1+2, 1+2)) 

00  139  J=1,NFS 
PFP  (I  ,  J)*PFP  (I,  J)  *VF  ACTOR 
PFP(J  ,  1)  *PFP  ( J,  I)  *VF  ACTOR 
PFP(I*2,J)*PFP(I+2,  J)*AFACTOR 
PFPCJ ,I*2)*PFP( J, 1+2)  *AF  ACTOR 

139  CONTINUE 

140  CONTINUE 
HANIN0*1 
GO  TO  121 

112  CONTINUE 
MANINO*0 

TRXXT* «8* TRXX TO 2* TRXXT 
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CALL  MAOO(XFP,DX,XFM,NFS,ONE,ONEI 
C  XFP=X(TI+)a  °<TI+)  *  R-l  *  HT  *  (Z  -  SHALL  H)  ♦  X 

IFCIFRINT.NE.l)  GO  TO  1006 
IF(TIME.LT.3.8>  GO  TO  1006 

PRINT  *,**XPEA<a”,XPEAK,"  YPEAK*",  YPEA  K,"  AR*",AR, 

♦"  SIGVF=“,SIGVF,"  TRXXTs**»  TRXXr 
PRINT  *,  -  “ 

PRINT  *,**XPP" 

CALL  MCUT(XFP,NFS,ONE> 

10  06  CONTINUE 
C 

C  ACQUISITION  OR  ESTIMATION  OF  QFO 

C 

IF<SIGF10.GT.0.)  GO  TO  74 
CALL  SHIFTA (QFO  »QF01 , NFS 2, NFS 2) 

IF (TI HE-«  19)  72,72,57 
57  IF(TIME-.5>  6i,8r,80 
61  IF  (TP  X XT -1000 • )  72,77,77 
C 

C  ACQUISITION  SCHEDULE  CHANGE  0*  QFO 

C 


72 


VARYQ=VARYQ-SIGF1*.  064'  44h444 
QFO(i,l>  =  DT**5*VARYQ/2'. 

0F0(l,3)  =  0T**4*VARYQ/8. 

GF0(l,5>=0T**3*VAPYQ/6. 

QFO(2 , 2)*QFD  (1,1) 

OFO  C2 ,4)  *  OFD  ( 1 ,  3) 

OFO  (2 ,6)  sQFD  (1,5) 

QF0(3 ,l)sQFO (1, 3) 

QFO(3,3>*OT**3*VAR'Q/3. 

QF0(3,5)  =  DT**2*VARY*>/2. 

QFO  (4, 2)*  QFO  (7,  It 
QFO  (4  ,4)* QFO  (3,  3) 

QFO  (4  ,6) a  QFO  (3, 5) 

UFO  (5 , 1)  =  QF0(1,  5) 

QFOC5 , 3)  sQFD  (3,5) 

QFO(5,5)*VARYQ*QT 
QFO (6 , 2)  sQFD  (1,5) 

OF  0(6, 4)* QFO (3, 5) 

QFO (6 ,6)  -  QFO  (5, 5) 

QFO  (7,7)  *  (SIGF2**2)*(1.-EXP(2.*0EIT/FTAU2)) 
QFO  ( 8,8) *OFO  (7,7) 

GO  TO  74 


C 

C  ESTIMATION  OF  QFO 

C 

77  PRINT  *,**  - 

PRINT  ADAPTION  STARTED  AT  ~,riNE,"  SEC," 

80  CALL  MA30(QF0,0X0XT,PFP,NFS,NFS,3NE) 

CALL  MAOO(QFO,QFO,PPFP,NFS,NFS,-t) 

C  QFO*  OXOXT  ♦  P (TI+)  -  PHI  •  PCTI-l ) ♦  •  PHIT 

C 

C  60UN0ING  QFO 

c 

00  84  1*1, NFS 
QFACT0R*1. 

IF  (QFO  (I,  I)  *GT •  1«  )  GO  TO  86 
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QF ACT  OR* 3* 

0F0a,I)*.l 
60  TO  65 
86  CONTI  NL'E 

IF (S0RT(QF0 (1,1) ) *GT  *QF0MAX( I)  )  3FACT0R=0F0MAX  (D/SQRT  (QFOCI, I) ) 
IF(QFACTOR.3E.i*>  30  TO  84 
0F0  (I  ,  I) =QFACT0R**2*QF0 (I ,1) 

85  CONTINUE 

00  82  J*1 «NFS 
IFCI.EQ.J)  GO  TO  82 
QFOCI,  J)*QFD(I,J)*')F  ACTOR 

62  CONTINUE 

00  83  K*1  f  NFS 

IF (I«  EQ.K)  60  TO  83 

OFD  CK  »I)*QFO  CKf  I)  *QF  ACTOR 

63  CONTINUE 

64  CONTINUE  • 

00  73  I*lf NFS 
00  73  J*1 »  NFS 

73  OFO(I  *  J)  *  *2*QF0(I,  J) ♦*8*0F01 (I, J) 

74  CONTINUE 

WRITE  DATA  TO  FILE  TAPES 

SAVE(l)  *  XS(1) 

SA VE(  2)  *  UT(ltl) 

SA VE<  3)  *  XS  (2) 

SAVEC4)  *  UT (2, 1) 

SAVE<5)  *  XFP(l) 

SAVE! 6)  *  XFP<3> 

SA VE<7)  *  XPP(2) 

SA VE( 6)  *  XFP (4) 

SAVE (9)  *  XCENTR 
SAVE(IC)  *  YCENTR 
SAVE(ll)  *  PFP(lfl) 

SA VE( 12)  *  PF3 ( 3,  3) 

SA VE( 13)  *  PFP ( 2» 2) 

SA VE( 14)  *  PFP  <4, 4) 

WRITE  (6)  SAVE 
60  TO  50 

101  PRINT  ♦,  "LOST  TRACK,  X  CHANNEL,  HEAS  CALLEO  “,IHEAS,"  TIMES* 
♦RUN  ",L 

60  TO  105 

102  PRINT  ♦,  "LOST  TRACK,  Y  CHANNEL,  MEAS  CALLEO  ",IMEAS,"  TIMES. 
♦RUN  ",L 

105  00  110  J*l,14  v 

110  SA VE( J)*0 • 

106  WRITE (6)  SAVE 
TIME*TIME+OT 

IF(TIHE.GT.TFINAL)  60  TO  99 
60  TO  106 
99  CONTINUE 

STOP  "FINISH" 

ENO 
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6  State  Filter 

PR06F AH  THESIS (lNPUT*/30, OUTPUT, r*PS5*INPUT,TAOE£sOUTPUT,TAPE8> 
COHMON/FLIR/  XFOV, YF0V,I HAX,NPIX,  SIGHS, SI GMF, SIGMA B, SI GFLR,RF 

*  ,  ASPRC  ,FI  MAX  ,'JT  (2,1)  ,IHEAS,AR ,XFN  (6)  , SIGVF ,PL2P (64  ,2)  ,  VMAX,SIGMF) 
*, RANGED, RANGE 

INTEGER  ONE 
REAL  I MAX 

OX  HEN  SION  Z(64),PHI(8,6),Q(3,3), 

♦  WORK (3, 3) ,W(3) ,TEMP(8,8) , TEMPI (f , 3)  ,  SAVE( 14) ,OFDl (6, 6) , 

*00  (8, 6) »PHIF (3,6) ,QFO (f , 6) ,SQQO (8 ,8) , XS (E) ,H(£- ,5) ,PFP(6,6) , 
*XFP(6) ,PFM(6,5) ,HF  (S4)  » EXTRA (6,6) ,RIH (64 ,1) , CXOXT (6,6) ,PPFP(6,6) , 
*HT<6,64)  ,WKAREA(50,5  3) , C (5) , R (64,  64)  ,RN(£*,64)  ,80(8,2)  ,PHIFT(6,6), 
•OF DMA  X  (6) 

WRITE  (6,1) 

FORMAT  (1H1) 


READ  AND  ECHO  OATA 
RE AO  ♦,  SIGS1 

PRINT  *,  “RMS  OYNAMICS  FOR  TRUTH  MODEL,  SIGS1  *  ",SIGS1 

REAO  *,  SIGMAS 

PRINT  *,  “RMS  TRUTH  MODEL  PACKGR3UN3  NOISE,  SIGMAS  *  **, SIGMAS 
REAO  *,  SIGFLR 

PRINT  *,  "RMS  TRUTH  MODEL  FLIR  NDISE,  SIGFLR  *  ", SIGFLR 

READ  *,  IMAX 

PRINT  ♦,  "TRUTH  HOOEL  HAX  INTENSITY,  IMAX  *  ",IMAX 

REAO  *,  SIGAT 

PRINT  *,  "RMS  ATMOSPHERICS  FOR  T*UTH  MODEL,  SIGAT  *  ", SIGAT 
REAO  *,  NRUN 

PRINT  *,  "NUMSER  CF  MONTE  CARLO  ?UNS,  NRUN  *  ",NRUN 

REAO  ♦,  T FINAL 

PRINT  ♦,  "FINAL  TIME,  TFINAL  *  ",TFINAL 

REAO  *, SIGHS 

PRINT  *,  "INITIAL  RHS  TRUTH  MOOE.  SIGMA  PERVEL,  SIGHS  »  ", SIGHS 
REAO  ♦,  ASPRO 

PRINT  ♦,  "TARGET  ASPECT  RATIO,  ASPRO  *  ", ASPRO 

REAO  *,XQ 

PRINT  *,  "INITIAL  X  POSITION,  XC  *  ",X3 

REAO  *  ,Y0 

PRINT  *,  "INITIAL  Y  POSITION,  YO  *  ", YQ 

REAO  *,ZG 

PRINT  *,  "INITIAL  Z  POSITION,  ZC  *  ",Z0 

REAO  *,XOOTO 

PRINT  *,  "INITIAL  VELOCITY,  X  OIRECTION  ,  XOOT b  *  ",XOOTO 

REAO  *,VOOTO 

PRINT  *,  "INITIAL  VELOCITY,  Y  OIRECTION  ,  YOOTO  *  ",YOOTO 

REAO  *,ZOOTO 

PRINT  *,  "INITIAL  VELOCITY,  Z  OIRECTION  ,  ZOOTO  *  ",ZOOTO 

REAO  *,ISPTL 

PRINT  *,  "SPATIAL  NOISFI  1  YES,  O-NO  ISPTL  *  ",ISPTL 

IF  (IS  PTL*  NE«  1)  GO  TO  2 
REAO  *,C(1),C(2),C(3) ,C(4),C(5) 

PRINT  ♦,  "SPATIAL  NOISE  CORRELATION  COEFFICIENTS!  " 

PRINT  *,C(1)  ,C(2)  ,C(3)  ,C(4),C(5) 

2  CONTINUE 

PRINT  *,"  " 

REAO  *,  SIGMFO 

PRINT  *,  "INITIAL  FILTER  SIGMA  VELOCITY,  SIGHFO  «  ", SIGMFO 
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RE  AO  *,  ARO 

PRINT  *,  “INITIAL  FILTER  ASPECT  9ATI0,  ARO  *  “,ARU 

REAO  *,  SI6F2 

PRINT  •,  “RMS  FILTER  ATMOSPHERIC  NOISE,  SIGF2  *  “,SI6F2 

REAO  *,  RF 

PRINT  *,  “FILTER  MEASUPEMENT  NOISE  RMS,  RF  *  ",RF 

REAO  *,  SIGF1Q 

PRINT  *,  “INITIAL  RMS  DYNAMICS  FDR  FILTER,  SIGF10  *  “,SIGF19 
REAO  *,  FIMAX3 

PRINT  *,  “FILTER  MAX  INTENSITY,  FIMAXO  *  “,FIHAX9 

PARAMETER  VALUES 

YFOV*8. 

NP I X*  8 
AT  A Ul*  14. 14 
ATAU2*659.5 
FTAU2*1./ATAU1 
OT  *  1./30. 

IREF*  1 

FI  MAX «ABS (FIMAXO) 

SIGVF*ABS  (SIGMF3) 

SIGFlsABS  (SIGF10) 

OF  DMA  X  (1)  *2. 

OFOMA  X (2) =OFDMAX ( 1) 

OF  DMA  X  ( 3 )  =40  0  . 

OFOMA  X (4) =OFDMAX (3) 

OFOMA  X (5)  *,5 
OFOMA X (6)  *.5 

INITIALIZE  TRUTH  MOCEL  VARIABLES 

CALL  RANS ET<75682) 

ONE  *  1 
NPS  *  8 
NMS  *  NPIX**2 
NFS*6 

NFS2*NFS*NFS 
NFSM2*NFS-2 
NIS  *  3 
RI*1./RF 
SN*SI GMA9/IMAX 
IFtSN.LE.O.)  3N*»  001 
SN*1./SN 

RANGE  C*SQRT(X0**2*Y!)**2*Z0**2> 

AGAIN  *  ,351096534  *  SIGAT 
IFILE*6 
OELT  *  -l.*DT 
00  5  1*1,8 
BO (I, 1)*0 • 

BO (I, 2)*Q , 

00  5  J*l,  8 
00  (I,  J)  *  0, 

SQQO<  I,J)  *  0, 

PHKI , J)  *  0, 

5  CONTINUE 

FACT*  (AGAIN** 2)  *( ATAU1**2)*( ATAU2**4) 

H<r 
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FACTi  *  ATAU1-ATAU2 
FACT2  *  ATAU1+ATAU2 
FACTS  *  2.*ATAU2 
XF0V*8. 

61  *  FACT/ (FACT 1**4) 

62  *  FACT /(FACT  1**3) 

63  *  FACT /(FACT 1**2) 

R1  *  1.-  EXP(2.*-ATAU1'0ELT) 

R2  *  1.-  EXP (FACT2* CELT) 

R3  *  1.-  EXP ( 2 • *  AT  A  U2*  0 ELT ) 

R4  *  QT*EXP(OctT*  FACT 2) 

R5  *  DT*EXP(2«*AT AU2*r>ELT) 

FILL  OUT  TRUTH  MOOEL  PHI  MATRIX* 

SEE  MERCIER'S  THESIS  FOR  DERIVATION 

PH  1(1,1)=  1. 

PHI  (2  »  2)  *  PHI(1,1) 

PHI  (3  *  3)  =  EXP(ATAUl*OFLT) 

PH  I  (4 , 4)  *  EXP(  AT«U2*DEIT) 

PHI  (4, 5)  «  0ELT*PHI(4,/J 
PHI  (5  ,5)  a  PHI  (4,4) 

PHI (6 ,6)  *  PHI (3, 3) 

PHI (7 ,7)  a  PHI (4,  M 
PHI  (7 , 8)  a  PHI  (4,  5) 

PHI(6,e>  a  PHI  (5,5) 

WRITE (6,11) 

11  FORMAT (///2X, "THE  TRUTH  MOOEL  STATE  TRANSITION  MATRIX  IS*"/) 
CALL  MOUT (PHI,NPS,NPS) 

FILL  OUT  DISCRETE  INPUT  MATRK 

BO (1,1) =  OT 
BO  (2,  2)a  OT 
WRITE (6,15) 

15  FORMA T(///2X, "THE  TRUTH  MOOEL  IN»UT  MATRIX  ISt"/) 

CALL  MOUT  (90,NPS,2) 

FILL  THE  QO  MATRIX  WITH  VALUES  USING  EXACT  INTE6RATI0N 
SEE  MERCIER'S  THESIS  FOR  DERIVATION 

00(1,  l)a  SI6S1 
00(2,2)  a  Q0(1,1) 

00(3,3)  a  (61*R1) /(2.*ATAU1) 

00(3,4)  a  R2*(62/FACT2**2-Gl/FACr2)-R4*62/FACT2 
00(3,5)  *  G2*R2/FACT  2 
00(4,3)  a  00(3,4) 

00(4,4)  a  R3* (Gl/FACT3-2**S2/FACr  3** 2 *2*  *63/FACT3**3) • 

*  R5*( 62/ATAU2*G3*  OT/FACT  3-2»*G3/rACT3**2) 

00(4,5)  «  R3* (63/ FACTS** 2*62/ FACT  3)  5*63/ FACTS 

00(5,3)  a  00(3,5) 

00(5,4)  a  Q0(4,5) 

00(5,5)  *  R3*53/FACT3 

00  2C  1*3,5 
00  20  J*3,5 
00(1*3, J*3>  *  00(1,  J) 

Q( I“2, J“2)  a  QO  (I , J) 
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20  CONTINUE  -•>  * 

WHITE  <6, 30) 

30  FORMA  T<///2X, "THE  TRUT**  MOOEL  QO  MATRIX  ISl “/) 

CALL  MOOT (  QD,NPS,NDS) 

,  c 

(  C  TAKING  CHOLESKY  SQUARE  ROOT  0*  00 

C 

SQQOUyl)  *  SORT  (00(1,1) ) 

SQ00<2,2>  *  SOQO(lyl) 

CALL  CHOLESKd,  WORK,  NIS) 

00  33  1=1, NIS 

00  33  J*1,NIS 

SOQO(  I  *2,  Jf  2)  *  UORK(J,I) 

SQQ0(I+5yJ+5)  *  WORK(J,I) 

33  CONTINUE 
WRITE (6,35) 

35  FORMAT  (///2X, “THE  CHOLESKY  SQUARE  ROOT  OF  QO  IS*"/) 
CALL  MOUT (SQOO, NPSyMPS) 

IF  CISPTL.NE.l)  GO  TO  M  • 

C 

C  SET  UP  SPATIAL  NOISF  CORRELATION  COEFFICIENT  MATRIX 

C 

N*64 

M*8 

00  36  1*1, N 
R<I,I)*1. 

IF  <1.GE.64)  GO  TO  36 
RU,I  ♦1>*C<1) 

IF  (I.GE.63)  GO  TO  36 
R(Iyl42)sC(3) 

f  IF  (I  .GE.  57)  GO  TO  36 

^  R(I,l4fc)*C<4> 

P(I,I*7)«C(2) 

Rllyl 48)=C(1) 

Rf  I,l49)*C<2> 

RCI,l4iO)*Ct4> 

IF  (I.GE.49)  GO  TO  36 
R(Itl4l4)*C<5> 

R(  1 ,1  415  )  *C  (4) 

R(Iyl4l6)*C(3) 

R(Iyl4l7)*C(4) 

RCI,I418)»C(5) 

36  CONTINUE 

00  37  1*1, M 
R(8*1-7,8*I>*3.0 
R(8*I-7,6*I-1)*1.C 
Rt  841-6,  8*1)  *3.9 
IF  (I.GE.8)  GO  TO  37 
RC  8*1, 8*1 41)  *3.9 
R(fl*I,fl*l42)«J.3 
Rf  6*I*1»8*I41)*0.  0 
R(8*I -7, 8*147) *0.9 
RC  6*1  -7,8*l48)*0.a 
R< 8*1 *6, 8*I*8) *3. 0 
IF  (I.GE.7)  GO  TO  3T 
R<  8*1,8*149)  *0.0 

;  Rf  8*1  ,8*1 4101*0  «C 
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Rce*i-i»e*i+9)ai.o 
IF  (Z *G£«6)  GO  TO  17 
R(fe*I,e*I+i7)*3.D 
R(8*l ,8*1 *18) *0  ,C 
R(8*I-l,8*I*l7)*0.3 

37  CONTINUE 

00  38  1*1, N 
L»Itl 

00  38  J*L,N 
R(J,I)*R(I,J> 

38  CONTINUE 

00  39  1*1, N 
00  39  J*1,N 
RN (I, J)*S1GMAB*R(I,  I) 

39  CONTINUE 
WRITE  (6,44) 

44  FORMAT  (///2X, "64  X  64  FPATIAL  CORRELATION  MATRIX*"/) 
CALL  M WRITE  (RN, 64,64,84) 

C 

C  COMPUTE  THE  CHOLESKY  SQUARE  ROOT  OF  RN 

C 

CALL  CK0LESK(RN,R,64) 

WRITE  (6,48) 

48  FORMAT (///2X, “THE  CHOLFSKY  SQUARE  ROOT  OF  RN*"/) 

CALL  N WRITE  (R,  64,64,  60 
60  TO  42 

41  00  43  1*1,64 

43  R(I,1)*1. 

42  CONTINUE 
C 

C  SET  UP  FILTER  MATRICES, 

C 

00  51  1*1, NFS 

00  51  J*1,NFS 

PHIF(I,J) *0, 

51  0F0(I  ,  J)  *0, 

C 

C  FILL  OUT  FILTER  PHI  MATRIX 

C 

PHIF(l,l)*i. 

PHIF( 2,2) *PHIF (1, II 
PHIF(1,3)*0T 
PHIF( 2,4) *0T 
PH  IF  ( 3,3)  *PHIF  ( 1,  1) 

PH  IF  (  4 , 4 )  *PHI  F  ( 1 ,  1) 

PHIF(5,5)*  EXP( 0ELT/FTAU2) 

PH IF (6 ,6) *PHIF (5,5) 

00  5€  1*1, NFS 
00  5t  J*1 ,  NFS 
56  PHIFT  (I,J)*PHIF (J,I) 

C 

C  FILL  OUT  FILTER  DISCRETE  QFO  MATRIX 

C  FOR  START  OF  ACQUISITION  PHASE 

C 

0F0(l,l)*0T**3*SIGFl/3. 

OFO  (  2,2)  *OFD  (1,1) 

OFO  (1,3) *DT**2*SISFi/2, 
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0F0  <3,1)*QF0  (1,3)  *•'  * 

0F0  (2,4)  »QF0  (1,3) 

0F0  (4,2)  «0F0  (1,3) 

QFD  (3,3)*0T*SI6F1 
OFD  (4,4)*QF0  (3,3) 

QFO  (5,5)  *  (SIGF2**2)Mi.-EXP(2.*DELT/FTAU2)) 

QFD  (b,6)*QF0  (5,5) 

WRITE  (6,40) 

FORMA  T(///2X,  “THE  FILTER  STATE  TRANSITION  MATRIX  IS***/) 
CALL  MCUT(PHIF,NFS,NFS) 

WRITE  (6,45) 

45  FORMAT  (///2X, -THE  INITIAL  FILTER  QO  MATRIX  IS*-/) 

CALL  MOUT  (  QFO, NFS»NeS) 

PRINT  *,**  ** 

PRINT  - 


PRINT  *,* 


REGIN  THE  MO*TE  CARLO  SIMULATION 


00  99  L*i  ,NRUf( 

TIME  *  0. 

XCENT R*0. 

YC  ENT  R*0  • 

FI MIN  *0 • 

IN  EAS  *  G 

FI  MAX* ASS  (FIMAXQ) 

AR*ARO 

SIGVF*A9S  (SI  GIF  3) 

SIGF1*ABS  (SIGF10) 

VARYQsABS (SIGFlO) 

RESET  INITIAL  COMOITIONS  FOR  MEW  RUN 

00  46  1*1  ,NPS 
XSU)  *  J  . 

46  CONTINUE 

00  47  1*1, NFS 
XFP(I)  *  0* 

XFM(I)*0. 

00  47  J*1,NFS 
OF  0(1  ,J)*0. 

PFP(I,J>  *  0. 

47  CONTINUE 

PROVIDE  FILTER  WITH  INITIAL  VELOCITY  ANO  POSITION 

RHOR*  (X0**2»2<1**2) 

RANGE* (RhOR*V3**2) 

XFP(3) *(2Q*X00T  0-X9*  200TQ)/(RH0R*«  330  02) 

RHOR* SORT  (RHOR) 

XFP(4)*(RHOR*YOOTO«YO*(  (XO*XOOTOf  ?0*ZUOTCI  /RHOR)  )/ (RANGE*. 00002) 
FILL  IN  P*  AT  TINE  0 


PFP(1,1)*25. 

PFP(2 ,2)*PFP(1, 1) 
PFP(3,3)*2Q00. 

PFP  (4  ,4)  *PFP(3,  3) 


non  ooo  ooo  ooo  ooo 
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PFP(5,5>**2  \  ..  .. 

PFP(E,e)*PFP(5.5) 

FILL  IN  QFO  AT  TIME  0 

0F0(1  *1)  *DT**3*SIGFl/3« 

QFO  (2,2) *QFD  (i,l> 

OFO  ( i»3)  *0T**2*SIGe  i/2* 

0F0  (  3,1) *QF0  (lf  3) 

OFO  ( 2»4) *QF0  (1,3) 

0F0  (4,2)*QF0  (i,3) 

OFO  (3,3)*DT*SIGF1 
QFO  (4,4)  *QF0  (3,3) 

QFO  (5,5)  *  (SIGF2**2)M1*-EXP(2.*DELT/FTAU2)) 

QFO  (6»6) =QFO  (5,5) 

TIME  LOOP  STARTS  HEPE 

50  TIKE  *  TIME  ♦  OT 

IF(TIHE.GT.TFINAL)  GO  TO  99 

PERFORM  TRUTH  MOO  EL  SIMULATION 

VTIME*TIME-DT/2* 

XVLH*  *5D3  •♦TIHE+X3 
YVEH=-3uO •♦TIMELY  3 
ZVEH=Z0 
XOOT*-5flO. 

YOOT*-300. 

700T*0. 

RHOR* (XV£H**2+7  VEH**  2) 

RANGE* (RHOR*YVEH‘*2) 

UT (1, l)*(ZVEH*XD0T-XVEH*ZO0T) /(RNOR**03002> 

RHOR* SORT  (RHOR) 

UT (2, 1)* (RHOR* TOO T-Y  VEH* ( (XVEH*XDOT+7VEH*ZOOT) /RHOR) ) /(RANG 
GOC  02) 

RANGE  *SQRT  (RANGE) 

VH AX* SORT (XOOT**2*YOOTr*  2+ZOOT**2) /(RANGE* *60002) 
XC£NTR*XCENTR»UT ( 1,1  )*CD  (i,i) 

YCENT  fi*YCENTR*UT (2,1)* 60 (2, 2) 

CALL  NOISE  (NPSyW) 

CA LL  MMPY (TEMP,  SOOO, W, NPS,NPS ,ONt ) 

CALL  MMPY (TEMPI, PHI, >S,NPS,NPS,ONF) 

CALL  HA0D(XS, TEMP, TEMPI, NPS, ONE, 3NE) 

CALL  MHPY (TEMPI ,80,UT,NPS,2, ONE) 

CALL  MAOO (XS«X  S, TEMPI, NPS, ONE, ON*  ) 

FILTER  STATE  PROPAGATION* 

CALL  MMPY  (XFM,PHIF,XFP, NFS, NFS, ONE) 

FORM  CENTROIO  POSITION  AND  FILL  TRUTH  ARRAY 

XPEAK  *  XS(1)  ♦  XS(f)  4  XS(4)  -<FN(1) 

YPEAK  *  XS(2)  ♦  XS(6)  ♦  XS (7)  -<FM(2) 

IF (A8S (XPEAK) «GT*3.*ASPR0*SIGMS)  GO  TO  1C1 
IF (AGS (YPEAK) *GT • 3** ASFRO*SIGMS)  GO  TO  1P2 
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CALL  MEASCXPEAK,YPEAK, A,Z,R>  .•>  • 

C 

C  SEARCH  FOR  FIHAX  ANO  FIHIN 

C 

IF  ( FI  MAX i)  «GT  •  0  • )  GO  TO  59 
FIHAX*0. 

FIHIKxC. 

00  55  1*1, NMS 

IFCZ(I).GT. FIHAX)  FT  FAX  *  Z(I> 

55  IF(?(I)«LT • FIHIN)  FIHIN  *  ZCI) 

59  CONTINUE 
C 

C  FILTER  COVARIANCE  PROPAGATION 

C 

CALL  MMPY (EXTRA, PHIFfPFP, NFS, NFS,  NFS) 

CALL  MHPY (PPFP, EXTRA, PHIFT, NFS, N?S, NFS) 

C  PPFP*  PHI  *  P(TI-l) ♦  *  PHIT 

CALL  HA90(PFH,PPFP,<)FO,NFS,NFS,0<E) 

C  PFM  *  P(TI-)sPNI  *  P(TI*1)  ♦  *  PHIT  ♦  G  *  Q  *  GT 

C 

C  PERFORM  MEASUREMENT  UPDATE  FOR  THE  FILTER 

C  INVERSE  COVARIANCE  rORM 

C  FORM  FILTER  CENTROIP  POSITION  AND  FILL  OUT  NON  LINEAR 

C  SMALL  H.  CALCULATE  PARTIAL  SMSLL  H  PARTIAL  X 

C 

XPEAK  *  XFM(5) 

YPEAK  *  XFM (6) 

IF(CXFMC3)**2+XFM(t»)  **2>.EO.O.)  <FM(3)*.0C1 
CALL  MEASF(XPEAK, YPEAK, ONE, HF,H> 

00  61*  1*1,  NMS 
DO  6C  J*1,NFS 
HT(J,I>  *  H(I,J) 

60  CONTINUE 

CALL  MMPY(PFP,HT, H,NFS,NMS,NFS) 

C  PFP*  HT  *  H  (R-l  IS  SCALAR  A*0  MULTIPLIED  LATER) 

IOGT  *  0 

CALL  LINV2F(PFM,NFS, NFS, EXTRA, I93T,  UK AREA, IER) 

C  EXTRA*  P(TI-)-l 

00  65  1*1,  NFS 
00  65  J*1,NFS 

PFP(I  ,  J)*PFP(I,  J)*RI  ♦  EXTRA  (I,  J) 

65  CONTINUE 

C  PFP*  P(TI<-)-l*  HT  *  R-l  •  H  ♦  p<rr-)-i 

IOGT  *  0 

CALL  LINV2FC PFP , NFS, NFS, EXTRA ,105 TfHKAREA, IER) 

DO  70  1*1, NFS 
DO  70  J*1 ,  NFS 
PFP  (I ,  J)  *  £XTRA(  I,  J) 

70  CONTINUE 
C  PFP*P(TIO 

CALL  MAOO(RIH,Z,HF,NMS,ONE,NFS) 

C  RI HsRESIOUALS*  7-SMALL  H 

IFCSIGMFO.GT.O.)  GO  TO  62 
C 

C  COMPUTE  NEW  ESTIMATE  OF  SIGVF  ANO  AR 

C 


00  56  1*1, NMS 
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AR*  AR*.001*(RIH(I)*PL2P(I,1>>  •• 

SIGVF*SIGVF^«  30 1*  (RIH(I) *PL2P (I ,2) ) 

56  CONTINUE 
IF(AR.LE.l.)  AR*i« 

IF4SIGVF.LE.0.)  sigvf*sigmfo 

62  CONTINUE 

CALL  MMPY (EXTRA, MT,RIH, NFS, NMS,0ME> 

C  EXTRA*  HT  *  (Z-SMALl  H) 

CALL  HMPY(XFP,PFP,EXTR* ,NFS,NFS,3NE> 

C  XFps  P(TIO  *  MT  «  (Z  -  SMALL  HI 

00  7F  1*1, NFS 
XFP(1)*XFP(I)*RI 
75  CONTINUE 

C  XFPsOELTA  X  *  P(TI*>  *  HT  *  R-i  *  (Z  -  SMALL  H) 

00  78  1*1,  NFS 
CO  78  J*i ,NFS 
OXOXT (I,J)*XFPCI)*XPPC J) 

78  CONTINUE 

C  OXOXT*  (X (TI ♦)  -  X(TI->>  *  (X(TI+)  -  XCTI-DT 

TRXXT*C. 

00  79  1*1, NFS 

79  TRXX7 *TRXXT*DXDXT (I ,  I) 

CALL  MAOO (XFP, X FP  ,XeM,  l'FS»0NE,0N: ) 

XFPsXCTI*)*  P  (  TI  ♦ )  *  R-i  *  HT  *  (7  -  SMALL  HI  ♦  X(TI-) 

ACQUISITION  OR  ESTIMATION  OF  3F0 

IF(S1GF10.GT.3«>  GO  TO  74 
CALL  SHIFTA(QFD,QF01,NFS2,NFS2) 

IF (TI ME-* 19)  72,72,57 

57  IF  (TI HE-«  5)  61,8C,81 
IF  (TPXXT-lQuG#)  72,77,77 

ACQUISITION  SCHEDULE  CHANGE  OF  QFD 

VARY0*VARYQ-SIGF1*, 0641444444 
QF0C1  ,i)*0T**3*  VA°Y0/3  • 

0F0<2 , 2)  *QFO(l,  1) 

0FD(l,3)*DT**2*VARY0/2. 

QFO(3,l)  *QF0(1,  3) 

OF 0(2 ,4)sQF0(l»  3) 

OF 0(1  #2)*QF0(1,  3) 

QFO  (3 , 3)  *0T*  VARYC 
OF 0(4  ,4)*QF0(3,  3) 

QFO (5 ,5)*  (SIGF2**2)  *(1#-EXP(2.*0ELT/FTAU2) ) 

0F0<6 ,6)*QF0  (5,5) 

GO  TO  74 

ESTIMATION  OF  QFO 
77  PRINT  " 

PRINT  AOAPTION  STAPTEO  AT  ",TIM£,"  SE3*“ 

60  CALL  MAOO(QFO, OXOXT, PFP,NFS, NFS, 3NE) 

CALL  MAOO (QFO.OFO, PPFP,NFS,NFS,-t ) 

C  QFO*  OXOXT  ♦  P (T Z -  PHI  *  P(TI-l)^  *  PHIT 

00  73  1*1, NFS 
00  73  J*1,NFS 
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73  OF 0(1  »  J)a  .2*0FD(I, J )*.8*QF01 (I, J)  -  •• 

7k  CONTINUE 

a 

WRITE  DATA  TO  FILE  TAPE8 

SAVE(l)  =  XS(1) 

SAVE!  2)  *  UT  (1*1) 

SA VE( 3)  *  XS (2) 

SAVE (4)  =  UT  ( 2 , 1 ) 

SAVEC5)  =  XFP(i> 

SA VE( 6)  a  XFP(3) 

SA VE( 7)  *  XFP(2) 

SA VE(  8)  *  XFPU) 

SA VE(9)  *  XCENTR 
SAVE(IO)  «  YCEHTR 
SAVE(ll)  *  PFP(1,1) 

SAVE (12)  *  PFPC3, 3) 

SA VE( 13)  »  PFP (2,  2) 

SAVE(IA)  a  PF®(4,4) 

WRITE  (8)  SAVE 
GO  TO  50 

101  PRINT  *,  -LOST  TRACK,  X  CHANNEL,  HEAS  CALLED  ",IMEAS,-  TIMES 
♦RUN  **,L 

GO  TC  105 

102  PRINT  *,  -LOST  TRACK,  V  CHANNEL,  HEAS  CALLEO  ”,IMEAS,“  TIMES 
♦PUN  “,L 

105  DO  110  Jsl,  14 
110  SAVE(J)ao. 

106  WRITE  (8)  SAVE 
TI HEsTIME+OT 

IF(TIME.GT.TFINAL)  GO  TO  99 
GO  TO  106 
99  CONTINUE 

STOP  “FINISH- 
END 


4  State  Filter 

FRuG‘  A ►.  i  HESl$(IN°JT  =  /.-  0, OUTPUT  ,r*°E5slNPJT ,TAP£- *0UT?UT,TAPE8> 
tur:hLu/PcArt/  XP3/,f?0Vf  iMAXfrtPIX,  Ei3**S,SI3HF, SIGMAS,  SIGFLR.RF 
u»F*M.»X,  J.  (2>i)  ,  it' t AS 
IN  T  tt  E  F  v*  N  £ 

REAL  lf*MX 

CXr  t?  SxU.f  Ttt,a>,PHl  (72,72)  ,0(3,3)  ,RI  (0,0  ,T  A(A,t)  ,W&:,(oA,  1), 
4HvjrK(o,3i  fVKw).Ti^3(T2,7  2)flE.rtPl(?2,?2),SArft(i  )  ,W?2(72,i) , 

•CO  (.‘I  ,7  2)  ,»Hj.F(>)  ,ine  )  ,S0Q0(/2,7  2)  ,  XS  <  7  2  )  tH(  ,A)  ,FFPU,*> 
*,XFP(»  ),  PF.1  (-.,•»),  X r  M  (•• )  ,HF(bO  ,£<TRa(4,A>  ,  RlHU*,!)  ,rfO(o4,l)  , 

*K»  U,0.  ),W<Ar\£M  (3u,1j)  ,C  (a)  ,K<Sh,  £-)  , RN ( c i  ,bA ) , 30 (72, 2) 

KtA U  *  f  OTA'J 

PRINT  ’  ,  “iKUi’M  MOO£L  DYNAMICS  CD  RRELm  TI  CM  TIME,  DTAU  *  ",OVAU 
p£m0  *  ,  jIGSl 

Pi\ I NT  *,  “KrtS  OYnAMICS  FOR  itsUTH  MODEL,  SIGSl  *  ",SIGSi 

REaC  «,  SIGMA 3 

PRINT  ♦,  "<Ma  TRUTH  MOCEL  6ACKGRDUO  NOISE,  STGMA6  *  ",SI3MA9 
k£fQ  SlGFL* 

PRINT  -,  "r'U  TfUTH  HOTEL  FLIP.  NDI3I,  SiGFLR  =  ", SIGFlR 

nEaC  *,  *MAX 

PRIM  ♦,  "TVJTH  MODEL  t-AX  INTENSITY,  1MAX  s  ",IM*X 

.'.’cao  *  f  sn 

PUM  -,  "KUTH  ‘iOOEL  SIGNAL  Tu  <OLSE  RAIIU,  SN  s  ",SN 

Pt<4j  x  f  oOSA 

Prvj-NT  4,  "k/.G  jY;»£  MICS/aT  MOSPHERT  C  TCISE  RATIO,  SOSA  =  ",SDSA 

P.Cp  J  *  ,  NiVUN 

FxiNT  -,  **.<Un  )ER  Or  MONlc  CAxLO  tUNS,  NRUN  *  ** , NRUN 

kEau  *,  TFi.MWi. 

PR*M  *,  "FiNAw  TIME,  TFINmL  =  ",TFiNAL 

REAL*  *,  alGFc 

p-cam  *f  *v.ms  filter  atmospheric  noise,  sigf2  =  ",sigf2 

KtfcO  4  ,S±GM5 

Pk I wT  *,  “lNlfiA-  RMS  "RUTH  MCJE.  SIGMA  FExVEL,  LIGMS  =  M, SIGHS 
FEAO  •“■ ,  SlGMF 

P.wNl  »,  'Vis  FIL'ER  BEAMWIDTH,  SIGMF  =  ",SIGMF 

READ  ’ ,  *C'J'0 

PRINT  -,  "COrOELATION.  1  YES,  D  XO,  ICOrtR  =  ",*COkk 

>.EpC  *,  i»F 

PRINT  ♦,  **FxLTEp  MEASUT  EM£NT  NUI5E  RMS,  RF  =  ",RF 

*EmJ  *,  SlGFi 

Ft*,,*  »,  "*MS  DYNAMICS  FG*  FlLTE? ,  SIGFl  *  ",SI3F1 

F£aC  FTAU1 

P.<*i,T  *,  "FILTER  OfMAMICS  COr.REUTIDS  Tlhi,  FTAU1  =  ",FT*U1 
r £  AO  « ,  FlMAX 

print  *,  "Filter  «ax  intensity,  fimax  *  ",fimax 

REmu  *i  a$P<u 

PRInT  4,  "TARGET  A3®£CT  RATIO,  ASPRO  *  ",ASPRO 

k£mu  4t  YXw 

PRINT  *,  "INITIAL  VELOCITY  In  X  DIRiGTiON,  VXw  *  ",VXm 

*  , VY* 

PRINT  4,  "INITIAL  YELOTITY  in  Y  DIrEGTION,  VYC  *  ",VY5 

IF(VXc4V2>VYu*v2.£C.ti.)  V  XL  *  •  Out 
KE Ab  *$  DTlM 

PRINT  *,  "TRUTH  MODEL  NGiSt  COr»R:LATlON  TIME,  OTIM  *  ",OTIM 

KEAU  ♦flaPTL 

print  ♦,  "spatial  noise*  i-yes,  i-nd  •  isptl  «  ",isptl 

KeAO  *,C(i),C(2),C( 1) ,C (f),C(G) 

F*xNT  "SPATaa.  NOISE  CORrElhTIO.i  coefficients*  " 


PMMT  *y0  (1)  .3(  2)  ( 1)  ,C(t)  ,C(&> 

XFcVsfe. 

YFu V*fc  • 

AT  kU1  =  1h  .  1» 

A  f  »Ui.  aci  »  .’i 
F i  *-U«.=1»/AIhU1 

or  =  i*/5 j • 

I-UF=1 

INITIALIZE  T^UTH  MODEL  VMUA'T.ES 

CALL  RAMJET (123«?) 

U-\! 

1Y  =  L 

XM=u. I 

YMsiw^  L 

CUE  =  1 

NPb  -  72 

ini  o  =  UPi  X,r  *  2 

NF  b-* 

r.lL  =  Z 

Fxsl« /hF 

AG„Ii,  =  «  3C11  *  SIGSI/SOSa 

i  F  i  Lt  -  e- 
OcLl  =  -1  .-QT 
Uf  (it  1)  =  .  . 

L'  l  ( E  t  i ) = • 

00  j/  1  *  i  «  is 
60  (1,  1)  =  .  . 

BO  • 

OJ  i  J=l,w 
03  (if  0)  -  J. 

SQQ0( 1 *  J)  =  b* 

FHiCI  ,J)  s  w. 

CONTI  NUE 

FACT-  (AGAiU**  2)  M  ATAU1*  *2)*  UTA«)2***> 

FUTl  =  AT  AUl-ATA IJ2 

FacT*.  *  *»T*Jl*ATiU2 

FAlTZ  =  2»-*ATVJ2 

Gl  a  FACT/(FAjrif--A) 

G2  =  F*C»  / (F43T  i*  “  3)  . 

G3  a  FACT /(FA JT i*  '2) 

M  *  1.-  fcXP  ( 2  •  ATAIJ1  OtLT) 

R2  a  1.-  tXP(FA3T2mDFl  T ) 

K  3  s  j..-  EXP(2.-  ATAU2*  OELT) 

P-*  *  0T»cXP(JEL7*  FACT2) 

a  CT’tXPtE.-ATAUP'OrLTJ 

K9*CVP(0LLT/UT*M) 

P.<iM  *,  ‘Vo  =  M  t  rs> 

PHI  (1  y  1)  a  1. 

PHi  (2  »  2)  a  PHKlyl) 

FHl(3y3)  a  EXP(Af AUi^OELT) 

PHidyO  a  £XPCATAlJ2n0FLT) 

PHiC.  ft  )  a  DELT*=>HI<4,'  > 

FHl(i,  ,S>)  *  PHI  (  »,t  ) 

PHiU  y  L )  *  PHI  (3*3) 
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PHi(«|7)  *  PHI('4,h) 

PH*  (/  to)  =  PHI  (*♦,  F) 

PHI  (f  ,c)  =  pHX  (  3,  f  ) 

CJ  9  1  =  3,7 2 
i  FH1  (1  1 2)  *l\*j 

WtITt  It, 11) 

11  FOFfiM  <2X,“T^£  TRUTH  MODEL  STATE  TRhSSITIOn  MATRIX  IS*"//) 
F  win  -v*  “ 

63(1,1):  Oi 
EsU(2,  2)  =  OT 
Wwlf  (t,l!5) 

IS  F0uK*'M2X,"lHE  TRUTH  MODEL  INPUT  MATRIX  1  if//) 

CALL  MCUT (3U,NP3,2) 

FILL  T He  03  MATRIX  PITH  VAlUES  USING  EXACT  INTEGRATION 

CJ(1,1)  =  aiSSl 
Cj(2,2)  =  OU(i,l) 

C"3  (3,3)  =  (Gl*M) /(  2.*TTAU1) 

CQ(3,U  =  R2*  (G2/CAC7  2J  ♦  2-Gl/FACT  2)  -  R4*G2/FaCT  2 
CO  (3,1)  =  G«i*R2/ecCT2 
CO  (•»,  3)  =  00(3, V) 

00  (•. ,0  =  :v3MGj./F4CT3-2.*32/FACT3**2  +  2.*33/FACT3**3)  - 
*  RI-CGc/aI  AU2*G3r  Dr/FArT3-2.*G3/?6CT3*-2) 

C3(i,0)  =  .so*  (G3/FaCT3' *2-G2/FACf  3) -Rv- G3/FALT2 
°0(y,3)  =  03(3, j*> 

0)(C,H)  =  00(4,3) 

00  (6, i)  =  i^3* 3 3/FA C" 3 
00  2t  1=3, r 
DO  2l  J=3,*j 
CO (!♦ 3, J*3)  =  00(1,  J) 

0(1-2 ,J-2)  = 

2*3  CL  Mi  NOE 

73LUGS=S*GHA9Ml.~Vj**2.) 

DO  21  1=6,72 
1  00  l  I)=7D0UGi 

WrUTF  (c,3o) 

3j  FORM- T  (2X,"f HE  7\UTH  MCOEL  QO  MATRIX  IS*"//) 

6000(1,1)  =  S  3\T  (.03  (1,1)) 

SQCu<  2,2)  =  6000(1,1) 

CALL  CH0i.tS<n,N3F<,Nir) 

DO  33  1=1, MIS 
DO  33  J=x,NxS 
C0Qi)(  i +2*  J+2)  =  M  OR<  ( J,  1 ) 

SOGD ( I 45, J+5)  =  WGR<( J,l ) 

3  3  CONTINUE 

V=60hl  (ZUOU5S) 

DO  3*>  1*5,/  2 
4  6000(1, I)=Y 

WR*U  (t,35) 

55  Fjr.Mf  T  (2X  ,”f  ME  CHO.ESKY  SQUARE  ROOT  OF  00  IS*"//) 

V=60K'l  (S*GrfA'T) 

IF  (ISFTL.NE.t)  30  TO  » 1 

set  up  spatial  (oise  cortelation  chfficient  matrix 

N*o*» 
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M*t 

to  2«  x=i 

r(I ,2 )*i. 

IF  (<  iGwii  >)  SO  i  0  3c 

Mi,:  ♦!» -C < x) 

IF  (1.Gs:.*j3)  SO  ro  36 
rA  In  +  Z)*C(5) 

IF  S  J  fo  36 

M  i*I  ♦t)*C(*.) 

MI, 2  */)sC(Z) 

F(i,l  ♦c)3Ul) 
k(i 

»*  ♦9>=C<2) 

Mi,i+iL>  =C(*r> 

IF  (i«Gc*4i)  SO  TO  30 
M 1  »i  *m)  =  0(1  ) 
r. (  a,:  ♦!*»>  *D(-.) 
k  (  I  ,i  *1©)  -y  (  3) 

^  <  x  ,1  +i/)s«(‘  ) 

?.(i,I+l6)  *C(i  ) 

76  CJN'IiNUc 

uO  27  i*i»*i 
Mu*!  -7  tv^i)  =  j»  J 
F  (  u’  I  T(  i«  S 

R C  -t, c*i ) *i«  0 
IF  (1  .GE.d)  GO  TO  3* 

Mt*I  ,t**+i)*J.  • 

MfJ  J 

M  t**  =  “>•  i 

k(c4l-?,o*I  +  ?)* J.  - 

*.  ,  j *i*o)  =  1*  If 

M  t*i -t, a*I*t ) *  0«  * 

IF  (i.Gi.?>  GO  TO  3' 

F(w#l  ,c*i*a)*J«  J 

!%(  L*  1  ,  it*  «  *1« )  *«  •  I 

Me*!  *  j.  u 

IF  (...Get?)  SO  TO  3 7 

r. (  L*I  ,  c*i  ♦!/ )  *v  •  j 

M«*i ,w*i+la) *. •» 

m  *«.  • . 

37  COhlJMJi 

00  3v  isl  |N 
L*2  *1 

CO  It  J*L,N 
M  Jm  )  *x(  I,  J) 

3o  COhTlNUZ 

DO  3;  1*1,. 4 

OO  ih  J*  i  ,  ■> 

>..<  (I,  j)*y*x(i,J) 

33  CONTI  nUc 

FRIU7  '  ,  ■*  •• 

FtiMl  <•,-«  3/  6-,  30VAFIANCE  HmTUXI" 
Fm.41  -,**  ** 

PmNI  ” 

COMPUTE  The  CHOwiSKI  SUCARE  ROUOT  OF  RN 
CALL  CHJLEiK  (  <N,X,u  ♦) 
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GJ  TO  *.2  •* 

«*  1  CO  -,3  xsitu41 

*3  »•.(!,:  j*y 

s  2  CoN'i  x  NUE 

Phi  in  *»'*  " 

PrwrT.  '  OHO.ESKI  COUAkE  ROOT  OP  fcll*" 

PUw’  V  " 

fr.*M  M 

S;7  UP  FiLf£^  TIME  INVARIANT  iATUCES. 

PHiFll)  =  EXP10£LT/FTAU1) 

PH*F ( 2)  =  PH.F(l) 

PHxF( 3 )  =  £XP(u£.T/F7AU2) 

PH*F(  *»)  =  PHir  (3) 

CP  Ci  (i  )  *  (JIGPl*1  2)*a.-£XP(2.*0iLT/FTAUi)) 

PFb(2)  =  OFO(l) 

GFb  (3  )  =  CixGF2*'  2)  *<1.-EXP(2.'0£LT/FTmU2>  ) 

Cr  L  (‘  )  =  QFO(3) 

WkITL  (b,  1 3) 

FC»vH*.  1  t2X.“7H=.  FIIT”  STATE  TRANSITION  NATkiX  jIASONAL  Tt*NS  ARE") 
Caul  nCUT  <phx?,  j- >,  one) 

WiUTE 

«->  FOr.iip  1  (2<  »**i  Hi  r'  LTEr  CD  MAT  KiX  DIAGONAL  TukMS  4  F  E 1  ’*/  /  ) 

CAlc  t. CU  i(  OF  i,  NFS,  ONE) 

PUNl  >  3EG1N  THt  MO^TE  GAkLG  SIttULATIuN  ♦♦♦»*" 

OCl  bS  L*1,nkU  I 

1 1 1*,  t  =  V  • 

iFtfiC  =  l 
YNtW=U. 

XNtW*  L  • 

IUPUF  *•  =  j 
Xv  L  n1  Ksi« 

YClu'i  P=J. 

XiOOUI*L. 

YOEUh*v. 

KcStf  INITIAL  r,OlOITAONS  FO*  *EW  *UN 

GJ  LI  i»lf.iP5> 

X3  ( 1)  *  a. 

*■  c  Cu  i«Ti  f,U£ 

DO  -♦/  i*i»NFS 
XFFC)  =  0. 

Go  *♦?  j*i,(F$ 

PFPU  ,J)  = 

IFU.tO.J)  PFP(I,J)  *  .S 
7  CONTINUE 

Call  Sold )  , 

Call  MhPY  (WO,.i,Ni>,  3v,i  4,1) 

00  *»e  1*1,0-, 

Ad  Xi  (6*1)*EIGHA3*  WO(I) 

IP  (If  0K.< •  £i3«  1)  CALL  NEAF  (XW£M,  YN£W,NPS,RA,  R) 

»-  0  Txp.l  *  TiHi  ♦  JT 

IF  <  Ti  Pt,GT»7FI  nAi.  )  GO  TO  S3 
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IPL*G  *  i.FLAG*l 
C 

C  PfkFo^ri  V.-\ JTH  103EL  SIMULATE 

•« 

W 

AK®  (t  T  /2  • )  “Tl  ic 
UT 11, 1>=  VX_' £X9(A^) 

UT  U|J)a'/Y»*LXP((.V2i) 

Xw  tl»"i  r  s-XC  cCi  fvVUT  ( 1 1 1  )*  ’’3  ( 1|  1 ) 

YClN*»  **Yl  ai«Tk>UT(  2,  l)*-pDl2,2) 

CALL  NlH»E  (C|W) 

CALL  hCxbt  (t+iUjU 

CALL  Kf.PY  <W0,*,wC,,6A,b4,l> 

DO  i*.  1*1,  i 
52  H/LU,1)  =  W(1> 

DO  rc  1sa»*j4 

^3  Mi‘  L  t «  ♦& i x  )  ®Wt»  (a.  i  1  ) 

CALL  MfPY  ( T2MP,,nnn,w?2lNPS,MP3,3N£) 

CALL  Ht-Pf  (TLM3l|“HI|XS|NPS|NPS|0^r) 

XU*D*.  *\j(  (1*1) 

YO®l>*«  ♦UT  (  2»1) 

CALL  KA00(X3,T£  1= ,  Ti*;?! ,  NPS,  uhE,  3N£> 

C«lL  «HPY  (T£M3l|Hn,  JT|».Pi,2,0N£) 

CALL  >*«  )  )(Xi*XJ|T£^l,!*P5.|ONciONI) 

CALL  PiX;.0(X0,Y0|X"1,YM,iX,iY) 

CALL  blATMO(Xi,IX, IY ) 

Call  FIl..  (X3,IX,i  Y,  3:G*A9> 

XPt*<  =  XiU)  *  X3(3>  ♦  Xi(L)  -X3SUX 
YPtAK  =  XS(2)  «■  XS(j)  *  XS<71  -YJSUI 
IF  t*-  Bb(XPcAK)  «GT  »3»v  ASPrtO^SI&'lS)  50  Tu  It  1 
IF  ir  rOtYPEAO  .51  .3.  *AFPK0#  biGrtS)  53  TO  U*2 

Qu  >1  i = x  )Ot 

al  Zt^JsXSU+j) 

CAuL  KcAS  CXPlA<,YP£AK|NFS,2|K) 

C 

C  Ptr<F0‘<rt  00<*ElATION  T KA3K1NG 

U 

xF  ( xC  Ur  \ x)  GO  TO  5^«  3 
lUrliA 1 ®IOPOA1 +1 
C 

C  FILL  THE  TmRGET  AP.FAY 

C 

CALl  SHirTMZ,rA,>H3|l«F5> 

CALL  CD  A,NPi,NPt  ,  XC0<*,tYC3'i*,CXfCYI 

XClN'i  k*XULW*XCOR* 

YCtNT  ►*Yi,£M+Y0u3n 
IF(iLP0«'»  ..i£.I^Er)5  3  TC  5XC 
CALL  aHXFTAtZf AAf D<iS|NhS) 

XUtWsXCtNTR 

YMs.W*YCEnT* 

lUFOhl®* 

5  3  %*  v  CO  n  1 1  MJi 

c 

C  Pt kF0<vM  FIlTE?  STATF  AMD  CO^A?IAiD£  PA3PGGAT ION 

C 

XFMl)*j. 

XFP 12  )  *u  • 

00  ui  **1,MFS 
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C  TO  tot  J*x,NFS 

xf ci )  *  phif(  * >  - xp*(  l> 

PFr  (2  ,  J)  *  PFP(I, J) *PH:F<I)*PHiF{ J) 

IF  u.  to.  J)  PF1(x,J>  S  pfmc;,j>  ♦  :FD(I) 

u.  CuaIIkUE 

Ft  NFJr.H  rtEA.a'.KElENT  UFU4TE  FOR  THE  FILTER 

X3tAK  s  XFlCi)  ♦  XFH<?> 

YPtAK'  =  XFM(2)  ♦  Xri(U) 

CAlL  HiAxF  (XP£A<»  vd€ AK  ,  OmE,  HF  »  H) 

LJ  bi  I»L,.VlS 
CO  Ci  J*i  ,NF3 
HT  U,i>  s  H  (I »  J ) 

v  ■  continue 

Call  Mf.PY  (CFP,HT,  H,NFS,  NMi,  NFS) 

IDG1  =  (i  .• 

Call  LI  #72FCPFM»seS,NFEtEXTHA,ID3Tf-l<AKE&»xER) 
CO  ot  1*1 »NFS 
0<J  t:  J*x  t  NF  3 

PPM!  ,  J)sPrF<I,  J)  ’ll  ♦  tXTxA(I,J> 
l!  CONTINUE 
I'J  G'i  *  J 

CALL  L1R7 2F (PrP, nF3»  NFS ,  EXTkAjlQST  , XKAKEA, IEP) 
“JO  'ii  x*  1  ,  NFS 
CO  ir:  J*l»i«FS 
FF f-  U  ,  J)  a  EXTnA  <  x»  J) 

?  -  CjM2  NU£ 

CALL  NAOK  ixH,Z,HF,N.1P,ONc,NFS> 

CALL  rtf.PY  (EXTRA, HTf-*IH, NFS, NMS,OXE) 

Call  FMPY  (Xr'P, PF®, EXTR*  , NFS, NFS, 3NE) 

ou  Vt  i*:,nfs 

XFPC)  *  XFPiI)*RI  ♦  XFfi(I) 

7s>  CONTINUE 

XOSUl.aXOSUH+XFPd) 

YJLUf  *Y03UM«-X-P<2) 

SAVE(l)  *  X.»<1) 

SAVt(Z)  *  XS(S)  ♦  X3  (  4) 

SaV£<2>  *  X$(2) 

SAVL(f)  a  XS(J)  ♦  XS  (  7 ) 

SAV£(i>  *  XDIJM 
LAVlU)  a  XFP  (3) 

S AVt</>  *  YUjJd 
SAVL(c)  *  XF»U> 

SAVE<*s)  a  XCtUP 
3A Vt(  1 .)  *  YCEHTR 
Sa  VE(  11)  *  PF3  ( 1,  1) 

Sm VE(  1 2)  *  3F3(T,3) 

EAWtilJ)  *  PF®<2f  2) 

SAVE!  Is)  s  PF»(‘»,  x) 

WkIIl  (fe)  SAVc. 

GO  TC  51 

Ml  FKiNT  -,"LOST  TRACK,  X  CHANNEL,  RUN  “,L 
GO  TC  lMO 

t<5  2  PRINT  *  »**L0ST  TRhCK,  V  CHANNEL,  RUN  “*L 
lw»  00  llC  J*lili 
lit  SA  V&( J)*w  • 

(»' 
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I  Subroutines 

1  SUBROUTINE  ME4S  ( X»EAK,  YPEAK,  ISU3 ,Z,  R)  -  • 

COKMGN/FLIR/  XFDV ,YFOV, IMAX, NPIX,  SIGHS, SIGMF,SIGM.ABfSIGFLR, RF 
*,ASPP0,FIHAX,UT(2,1),IMEAS,AR,XF'U8)  ,SI6VF,PL2P(64»2)»VHAX,S1GMF 
♦,  RANGED,  RANGE,  SIGPVF,  C«?TH,SNTH 
REAL  IPAX 

(  DIMENSION  Z(8f 8)  i  R(6*»,6*>  , HI (6* ,1)  , HO <b4 ,1) 

ZHIN  =  0. 

SIGPVsSIGMS*RANGEO/RANGE 
PLVEL*SQR7  (UT < 1* 1) **2+UT  (2#1)  **2) 
tfl  SNlHsUT(2,i)/PLV?L 

CSTH=UT{1,1)/Pim 

SIGV* (!•♦ (ASPRO-i* ) *PLVEL/VMAX) *5IGPV 
I  «  <NFIX*ISU9> 

ID1V  a  ISU9**2 

15  XI NCR  *  XFOV/FLOAT(I) 

YINCF  *  YFOtf/FLOAT  (t) 

X  «  -li*XFOV/2.  ♦  XINCPV2. 

V  *  YFOV/2.  -  YINCR'2. 

Xfl  *  X 

29  CALL  N0ISE(6«,m) 

CALL  MHPY (H0«R*Wl«G4y6^  »1) 

MN=C 

00  2C  Ksi.NPIX 
00  13  J=1,NPIX 
2K  MN=MN*i 

TOTAL  =  C .  ' 

XN  =  X 
YN  *  Y 

00  1C  KslflSU? 

30  YCSTH*  (YN-YPEAK) '  CSTH 

YSNTH*  (YN-YPEAK)  *  SNTH 
DO  5  M  =1,ISU3 

ARGSP  VsYCSTH- (XN-XPEAK) *SNTH 
ARGSV* (XN-XPEAK) r  CSTH+YSNTH 

35  ARG«-  (  (ARGSV/SIGV) *  *2+ (ARGSPV/SISPV)  **  2)  *  •  5 

TOTAL  sTOTAL+£XP( ARG)*IMAX 
XN  *  X  ♦FLOAT (M)^XINCR 
b  CONTINUE 

YN  *  Y-  FLOAT <N)*YINC9 
40  XN  *  X 

10  CONTINUE 

U K,J) *TOTAL/FLOAT< IOIV) 

C 

C  AOO  BACKGROUND  MO  FLIR  NOISE  BOTH  ZERO  MEAN 

45  C 

IFCSIGFLR.EQ.O.)  GO  TO  30 
GA  USS  *  C • 

00  110  LL*lf  12 
GAUSS«GAUSS+RANF(0UH1> 

59  110  CONTINUE 

F*(GAUSS-6.)*SIGFLR 
ZCK,J)  *  ZCK.J)  ♦  F 
30  Z(K»J)*Z(K«J) ♦NO (MNfl) 

IF (Z( Kf J) »LT»ZMIM) ZNIN=Z  (K»  JJ 
5r  X  *  XN  ♦FLOAT (ISU9) *XINCR 

15  CONTINUE 

-  Y  *  Y  -  FLOAT (ISUP)*YINCR 
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xc 

CONTINUE 
IMEAS*IMEAS+1 
IF<7KlN.EQ.0.>RETURM 
00  25  1*1 »NPIX 
00  2S  J*1 »NPIX 
ZCI.J)  *  7<I. JJ-ZMIN+.l 
CONTINUE 
f  El  UP  N 
ENO 


subroutine  mmpy 


7«*/74  OPT  *1 


FTN  4.7+476 


SUBROUTINE  Mh»Y (C, A ,B ,K,H,N) 
OIMENSION  C<K»N)»A(<»M)fB(MfN) 

00  1  1*1. K 
00  1  J*1.N 
CU»J)*G. 

1  CONTINUE 
00  5  L*i.K 
DO  5  J=1.N 
00  6  1*1.  H 

Ca,J)  *  CCL.J)  ♦  (A(L,I)*B(If JJ» 
5  CONTINUE 
P.ETUP  N 
ENO 


SUBROUTINE  NOISE 
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FTN  4.7+476 


SUBROUTINE  NOISE <N,W> 
OIMENSION  W(N> 

00  15  J*1*N 
TOTAL *0. 

DO  5  1*1.12 

TOTAL  -TOTAL  ♦  RANF(OUN) 
5  CONTINUE 

W(J)  *  TOTAL  -  6. 

IS  CONTINUE 
RETURN 


SUBROUTINE  M&ASF(  XP£AK,  YPEAK,  ISU3 ,Z,  H)  ' 

COMMCN/FLIR/  XFOV,YFOV,lMAX,NPIX,  SIGHS, SISMFtSIG1AB,SIGFLi<fRF 
♦,ASPFO,FIKAX, UT(2,l),IMEAS,Aft,XFN(8)  , SIGVF,PL2P (64,2) , VMAX,SIGMF 
*,RAnGEO,RAnGE,SIGPVF,CSTH,$NTH 
DIMENSION  Z(e,8>,H(S4,f) 

REAL  I  MAX 
ZN1N  *  0. 

PLVEL  *  SORT  (XFM  (3)  **2+XFM(4)**2) 

SNTHa  XFM (A) /PLVEL 
CSTH=  XFM (3) /PLVEL 
SIGPVFaSIGVF/AR 
I  *  ( NP1X*ISU3) 

IOIV  «  ISU3**2 
XINCF  a  X FOV/FLOAT (IT 
YINCR  *  YFOV/FLOAT ( I) 

X  *  -l.*XF0V/2.  ♦  XINCF/2. 

Y  *  VFOV/2.  -  YINCR/2, 

X3  a  X 

DO  2i  K=1,NPIX 
NUM  a  K 
DO  lb  Jal.NPIX 
TOTAL  a 
SUM  1=0. 

SUM  2  s  Q  • 

XN  a  X 
YN  a  Y 

00  ll<  N=l» ISU3 
YCSTHa  (Yli-YPEAK) » CSTH 
YSNTHs(YN-YPEAK)*  SNTH 
DO  b  M  *1 , I SUB 

ARGSPV=YCSTH-(XN- XPEAK) ♦SNTH 
ARGSVa (XN-XPEAK)'  CSTH+YSNTH 

ARGa-  ((ARGSV/SIGVF) ‘*2+( ARGSPV/SIGPVF)*«2) *.S 
PART*  EXP (ARG) *  FI  MAX 
TOTAL  »  TOTAL* PART 

SUM1  a  SUM1  ♦  PAPTMARGSV*CSTH/SIGVF**2-6RGSPV  SNTH/SIGPVF**2) 
SUM 2  a  SUM2  ♦  PARTMARGSPV*CSTH/SIG»VF**2«-ARGSV*SNTH/SIGVFi-*2) 

XN  a  x  *FLOAT (M) *XIMCR 
5  CONTINUE 

YN  a  Y -FLOAT (N)  •  VI NC* 

XN  a  X 
IS  CONTINUE 

Z(KyJ)  *  TOTAL/FL OAT (IOIV) 

IF(SIGhFO.GT. J.)  GO  TO  16 

PL2P(K*(J-1)*3,1)*  PART* (-(ARGSP/ /SIGVF) **2*AR) 

PL2PC  K*(J-i) *8»2) a  PARY* ( (ARGSV**  2*4  RG3PV*' 2*AR**2) /SIGVF**3) 

16  CONTINUE 

IF (Z(KyJ) •LT.ZMIN)?MINaZ(K,J) 

H( NUM  $  1)  *  SUM1/FLOAT (IOIV) 

H( NUM  , 2)  a  SUM2/FL0AT ( IOIV) 

H(NUM,3)«0* 

H(  NUM  , 4) *  0* 

H(NUM,5)*Q. 

H(NUM,6)*0. 

H(NUH,7)*H(NUM,1) 

H(  NUM,  8)  *  H  (NUM  ,2) 

X  •  XN  ♦FL0AT(1SJB)*XINCR 
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j 

r 


NUM  a  NUH  ♦  NPIX 
15  CONTINUE 

V  *  Y  -  F10AT(1SU9>  *YIMCR 

x*xo 

20  CONTINUE 

IF ( ZM IN, EQ*0«) RETURN 
00  25  1*1 »NPIX 
00  2i  J*1 ,NPIX 
Zd,J)  *  ZCI. JI-ZMIN+.l 
25  CONTINUE 
RETURN 
ENO 


SUBROUTINE  HOUT  74/74  0PT=1 


FTN 


1  SUBROUTINE  MOJT  (6  ,N.R,  NC> 

01 HEKS10N  A  (NR,  NO 
00  It  1*1, NR 

WRITE  (6,5)  (A (I, J) » J*1»NC) 
5  5  FORMAT (2X,d(G13«7,3X)) 

10  CONTINUE 
RETURN 

,  ENO 


SUBROUTINE  MADO  74/74  0PT*1 


FTN 


1  SUBROUTINE  MAOO (C  »A ,B, J,K,IFtAG) 

DIMENSION  A(J,K>,9(  J,K>  ,C(J,K> 
IF(IFLAG.EQ.l)  GO  TO  6 
00  5  N*l,  J 
5  00  5  M*1,K 

C(K,M }  *  A (N»MJ  -  3 (N , M) 

5  CONTINUE 
RETUFN 

6  00  IP  N*1 ,  J 

13  00  1C  H*l, < 

C(N,M)  *  A(N»D  ♦  9  (N, H) 

10  CONTINUE 
RETURN 
ENO 


4»  ?  4-476 
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»s 


»<* 
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SUBROUTINE  CH3LESK(A,S,N) 

01  KEN'S  I  ON  A (1)  « S(  1) 

NO  lHs  N 
NOIMlsN+1 
TOL*l.E-e 
HR*0 

NN*N*  NOIM 

TOLl*  P  • 

00  1  1*1,NN,N3I11 
R=ABS  (Ad)) 

1  IF(R.GT.T0L1)T011=R 

T0L1  *  T0L1*1.£-12 
11*1 

00  5C  1*1 »N 

IH1  *  1-1 

00  5  JJ*1 »NN, SOI) 

5  S(JJ)  *  u. 

10  =  lH-I  HI 

R* A (10) -DOT (111,5  (II) ,S  (II) ) 

IF(ABS(R)  ,LT.(TCL'A(I0)+T0L1)  )  G3  TO  53 
IF (R)  15,59,23 
15  MRs-1 

WRITE (6,1333) 

iODC  FOkMA T (37H  TRIEO  TO  FACTOR  AN  INDEFINATE  MATRIX  ) 
RETURN 

20  S(IO)  *  SQRT(rt) 

HR=MR  +1 

1FU.EG.N)  RETURN 

L*1I+N0I1 

00  21  JJ*L,NN,NOIH 

IJ*JJ*IM1 

25  S(IJ)  =  ( A  (I  J)  -DOT  (1 11 ,  S  ( II)  ,  S  ( JJ ) ) )  /  S(ID) 

SO  II*II ♦NOIM 
RETURN 
ENO 


PJNCTION  OOT  74/74  0?T*1  FTN  4.7+«*7b 

•c 

1  FUNCTION  OOT  (SR,  A,  4) 

DIMENSION  A (1)  ,9(1) 

OOT*C  • 

00  1  1*1,  NR 

5  1  COT  *  OOT  ♦Adi  *  3(1) 

RETURN 

END 


SUBROUTINE  SWIFT A  74/74  0?T*1 


PTS  4.7+h76 


SUBROUTINE  SHIFTA (A N) 
DIMENSION  A  ( N)  ,B(N) 

00  ICO  K=1,M 
lb  3  B(K)  *  A  (  K) 

RETURN 

END 


SUBROUTINE  RNOF  74/74  OPf*l 


FTN  4. 7+476 


SUBROUTINE  RNOF(A,<) 

C  SUBROUTINE  TO  ROUNOOFF  A, TO  NEAREST  INTEGER, K 

C 

K*IFIX(A) 

b*a-k 

IF  (B,  LT,C  ,5)  GO  TO  100 
K*K«-i 

103  RETURN 
ENO 


C 
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•J-^3!*TIN£  MW-IT£ 


74  >74 


0»7=1 


rTN  'o 


1  Sl"?KOU'.Ir,E  (lyNyMiIOI") 

CI"tN?10*i  X(  111*1,1) 

PRINT  *  y **  ** 

00  3‘:  0  1**1,  N 

5  WRIT?  (4(I'-y  JL)  y  J'-*l ,  M> 

340  FO  F M7  T  <1X,16(F7.«  ,1X> 1  , 

3C3  CONTINUE 

PRINT  \* 

END 


S'J'XR^UTINE  PILL  7-/74  0*»T=1 


PIN  4 


1 


3 


3 


3 


a 


<4/1 


SUHHC  UTINE  PILL  (  X"  ,  I V  ,  IV  ,SIG«A,,> 

c-r .■>■£*■  sion  vs (i)  ,w m 

P=FLf  AT  ( I  X) 

?=FLCA7  (IY) 

N*£*c  <R) 

MsAPS  (S) 

ip  ax)  2,ir,i5 

5  co  f.  J=1,N 

CALL  N0ISE(3tW) 

00  6  1*1.8 

6  XS  (6*  I  +  J)  sSI&  lAV  M(T) 

GO  TC  10 

15  no?  1*1, N 

CALL  NOISE (8  » W) 

L*«3-I 

00  7  J*l, 8 

7  XCI5'  J*L)*SIG1«VU(  J> 

10  CONTI NLE 

L*  C 

IF  (TY)  2%2t,2) 

20  P 0  2l  1*1, m 

L=9-l 

CALL  KCISEU,R> 

00  21  J*1 , 3 

21  XS  (*4  L+J)  *3IG*A3*M(  j) 

GO  TC'  23 

30  CO  Zf  1*1, m 

CALL  NCZSS(a«'«) 
ro  2*  J*  1 , 3 

26  XS  IS-  14J)  atZGIM*  M(  J) 

25  CONTINUE 

frETUf  N 
FMC 
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./ 


-TK  <*, 


( 

5 


t  3 


•> 


C 


3 


su?'.ruu><E  sT-\’,-u  (xf,tx,:y> 
oi'Ef  iior,  xsm 

F=i 

R  =  FLC.AT(IX) 

*=c-F  PS(  ’) 

:o  it  i=i,d 

Ha  L  - 1 

CO  It  J*l,< 

IF  <IX>  v  , «• , 6 
5  M=c-j 

X$  <N«.*)  =  XS<M*  t-I*) 

GO  •r0  13 

ft  Mr  J 

YS  <N40=XS(N*i  +  IX) 

13  CCM7HLE 

S=FLC A* (IY) 

r.=  jf -a  es  f  s  > 
k=i 

cc  21  J  =  1 ,  N 
TO  21  J=1,B 
Ir  CY!  It, If, 13  • 

1“  v=&-: 

Xs  fc~M 

xs (<*J)=XS(K+ J-l^’Y) 

00  TC  ?1 
lft  H=I 

Ksc1  “ 

XS(K*J)=XSK«-J*  VTy> 

21  CONTINUE 
FETUF  N 
END 


Si.n? 31^1  NE  PIX-0  7 U°7  =1 


c’TN  *.7*‘  'o 


3 


SUPrCUTiNE  PIX-T  (  3X  ,  T  Y  ,X«  ,  YN ,  TX  ,  I O 
C  0IVcN  3X  AN3  nY  t»  FIXE Lf/TI^F  UNIT 
C  ESTIMATE  X  AfiO  Y  VELOCITY  IN'  lNTE",:k  <G.  OF  PIXELS 
IX=C 
IY=C 

XM*y,’  *nx 

YHsY,-  4TY 
Tl  =  AF  s  (X  4 ) 

T2  =  AC  S (Y^) 

IF  (T  a  «G c •  ,3>  00  TO  i 
GO  TC  2 

1  CALL  RNOF  <XM,IX) 

2  XH  =  XK-FLOAT(IX) 

IF  ( TZ  .GE.  .fi)  00  TP  3  ' 

GO  TO  *• 

3  CALL  PNOF ( YY , t Y ) 

•*  VM*Y*-eLOAT(IY) 

FE7UF  K  70 

END 


A 


Appendix  G 


PLOT  Routines  Listing 

This  appendix  contains  the  FORTRAN  code  for  the  statis 
tical  computations  and  plotting  program  (PLOT)  used  for  all 
filter  performance  plots  in  this  research.  Depending  on 
the  variables  of  interest,  the  labels  can  be  (and  were) 
changed  and  appropriate  data  processed  to  yield  desired 
plots. 


Plot  Program 


"GRAM  PLOTTED  !u/lu  0‘>T=2 


FTN  -r.7+47 6 


F-CG*  f  P  PLOTTED  (OUTB'JT  ,TaPE8  ,  PLOT  ,TAF  E6=0UT  PUT) 

REAL  KalCO\R 

CO»*K  N  XM  (1»Z)  (1"2)  ,FVAR  <i.?2>  ,  XiP  VAR  1 1 1  2)  ,  XMMVAR  (152) 

CGMM(  N/*mTh1X/<  ALCOPR  (is  ,2u,  150,  STATS  U  ,3,1ft)  ,FRR0fc<2, 2,  151  )  , 
*  7£«r  (15  2) ,1 Y i£(5r2) , VPT  (S,L,  21) 

REWIND  8 
NF'=/ 

rrr=i./3o. 

Fc  Hi)  G  9 
N'FlLS*fc 
NTiHL  *15  * 

NR  UN“  2 1 

CALL  FPOCESS<  4F7,;4<UN,f  TIMP,D7,N-  ILS) 

CALL  CL'TPUTOiFs  ,*,rTMr,6,nT) 

CALL  G-APH<NTI*-,M<??JN) 

ST  C-P 
END 


riNE  OUTPUT  74/7«  UPT*2 


FTN  «,.7+s76 


SUBROUTINE  OUTPUT (UE$,NSTAT,NTIM£,UpIL£,rT) 

PEAL  KALCCaP 

CO  NHC  n/MATPIX/K  ALCORR(  14,20, 15u)  ,  STAT  S  (a  ,  3, If  C  )  ,  ERROR  ( 2, 2, 15l  )  , 

♦  ZERO  (152 )  ,TY'1E (ir>2)  ,  VPT  (4,s,21> 
no  5;  1*1, NES 

TIME*  I  • 

WRITE  (KFILE,1)  I 

1  fqk*7T  (//2X,m»***  »•'**••  THE  FOLLOW!  4G  J.S  OUTPU*  FOR  ERROR  STATE  ** 

*  II,  2X, 

W-UTc.  (NFI LS, 2) 

2  P0F  MA  T  ( ZX  ,**TiiE**»  C,X  ,**MCAN  ERROR**,  FX ,  "STD  OF  ERP0R**,5X ,  “FiLTER  EST 
«,“l*«fTE  OF  STT-/) 

00  2"  J«lvNTX1c 
TI*E  *  TIHE+OT 

WniTt  (KFILE,3)TIWE,  (STATS(I,K,J),w*l,N3TAT) 

3  FOrsHAT  (1X,FG.2,3(TX,E12.»)) 

25  CONTINUE 

*  v  CONTINUE 
FETUFN 


I 


i"»r?HTIut  PF.r  CeS  3 


./r 


TK 


C 

C 


c 


C 


SU?Wtnri£  PsJCSi.SC4PStNHUMf  NTI'T:  ,DT,  NFILE) 

'•cAu  KALCO?'3 

COf  rtC  N/M-.  T'tlX/K  U  CO?'5  ( 1  4, 2C,  1** .  ) ,  S  UT  S>  ( -  ,  3 , 1  sU )  ,  ERROR  (2,2,15.), 
-  *c.'  0(i>2),TY,'l*(i'?),VP|  (n,  ‘,2H 

fimSUN  ESU  I(')  >S*JMS''  <f,)  ,\MSU M‘  ) 

DJ  i=l  fNRU  * 

•30  5-  J=1,NU1E 

F.EAOtNFlLE)  KALI  *1*  (L ,  I ,  J)  ,  L*l,  1  ' > 

•  '  CONTINUE 

r= FLl A 7(  1RUN) 

Fm  =  FLO«r  (?»r.'JN-<  > 

03  1*1,  NT  Iis 

?Eru(i)  = 

TY.'E<i)  =  FLCATO'OT 
co  i  k=i,  : 

*r3(>“'(IO  s 
S'J'i^r  (<)  = 

1  CCITir.L" 

no  2  “’=1,'. 

K.isur  m  =  c . 

2  COMIl.Uc. 

CO  :  L  *1 , ^ 

TO  j  jsi.N.-'JM 

FRJ  =  <AwCO:  W(L»  J ,  T)  -  <ALCCRR(LK  , I) 
ti»U  *1  ( L  )  =  ES'JKL)  ♦  ER- 
fjrCC(L)  =  211 1o°  ( L '  ♦  r****2 

f  lilJf  <U  *  OTJH!.)  +  KALCG4K(L+trf  J,  I) 

•)  CONTINUE 

r0  1'  J=1,'^U4 

Fk.'  1  =  <ALCO-  <ti,  >,  T)  -  KALCOKKf;,  J,1 ) 
t»  i  2  =  KALCO*  <<3,  J,I)  -  KALCOR^lf  ,  J,I) 

FSUM<*  )  =  23U  1(0)  *  Eorl 

FCUMU)  *  ESU1CS)  +  F*F2 

^0^(1)  s  sUIS'H”)  ♦  rkRl*  *  2 

SUpSC(I-)  *  3Ui$0< •>)  ♦  FRR2<*2 

j.  com:  sue 

CO  2.  K*  1  ,*» 

E^csN  *  £S!JM (<)  /c 

F3T0FV  *  SCOT  (  <311*^0  IK ) -R*cM EAN»* ?) /S M) 

fsicev  =  so*-t  ( r.  ic,m  ( <)  /r.) 

S7a7f(K,ifi)  *  <-y.FAi 
fV  ATS  (K,  2  ,  t)  * 

$TaU  <K,  3,1)  *  *7T0-V 

z;  conisua 

00  21  K* l ,2 
E-.EaN  =  ESUM  (<♦  »)  /p 

ESTDtV  *  SOPTt  (TUISO  (<♦**>  -MEMrA4‘"2>  /  ?M) 

E«.fOf  ( K,  1,1)  =  Er*rVJ 
E'it  c^  of, 2,1)  «  srroEv 
21  CONTINUE 
<0  CONTINUE 
DO  o!  Lsl,* 

00  Cl  1*1,4 

ESLM<l)sKALC0V(l,l,I*33)~K4LC0?UL*<  ,1,1*30 
SUMSC  II)sEiUM(I)' '.? 
yi  CONTINUE 

VPT(l  ,1,1)  *  VPT  (L  ,  2 , 1 )  *  VRT(L,3,l)a  VPT(L, ~,l)a  (. 
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ItMJTlNt  PeCCcS3 


7t/7>  0Pr=2 


FTN 


CO  5:  1*2,21 

RsFLCATC  ) 

e'isn.cATci-i) 

DC  ii 

II  IMFaJ*33 

E\r  =**  LCORR  Cl ,  I  ,  I  It)  -KALCORR  (l>  *•  ,1,  ITIhE) 
FSUr'.C  J)*  ESU*(J)*E?.R 
SUhSrCJ)*  SUM30(J)+ERp'*2 
EH  EAI  =  £SUM(J)/-» 

ES*n?  V*  &QkT((:»U>(S"CJ)-R*EMEAI««>*2)/RH) 
VBTa,J,I)s  E3TOEV 
6*  roNTINUt 

comikue 
b=>  C3M7MJE 
RE TUP K 
5HD 


3"TI NF  PT«1  OPT*?  FT*  4.7*47o 

C 

SUAkL  UIInE  PTiIT  ( A  » REFLAG) 

P£  lL  F6LCC7R 

CO  KMC  A  Xu  (lr  2)  »  V*  p  ( 1*12)  ,  FV®R(lr  2>  ,  X1F  V  AR  (1*  2)  ,XH«VAR  (152) 

COf,MC N/M4T<IX/*Al  pO  V?  ( 14 , 21 ,  IS O »  STATS  <l.  ,3, 150  ,  ERROR ( 2, 2, 150 
*  TEKt  (152),TYHE(1*'»)  ,  VcT  (**,4 , 211 
DIMENSION  A(m,«U17> 

VAF  (lel)*FVAP(l  ;i)*XMPVAR(l51)*X-HVi;  (151)*'’, 

V4r»(l*.  2)  *FVAK  (1*3)  *XMPVAR(1J2)  *XiHVAF  (152)  *C. 

CALL  SCAL E ( XHP W. , 4 •  ^ * 3 i» 2 »1) 

XHFV*  fi  (1*>1)  *X  1(1,1)  *?EP0(1*  1)*XH1V1R(1.11) 

VHPVAP  (1^  2)  *X  1  ( l:  2)  *’EP0(1‘  2)*XH1V4R(15  2)*-XHMVAS  (152) 

CALL  VCR*iPM(TrM-,y-,l5:  ,  A, -1,15,4  ) 

CALL  VGRAPH(TYM",Xi3VAc , 15J , A, 2,IE , 1‘ ) 

CALL  VGRAPH(TYM?,XM1VAr,  15P , A  , 2,1£  ,  14) 

IF (IF  LAG*  E0«  p) GO  *0  qq 
CALL  SCALE(VAl,U.p,  132,1) 

VAR  (it  1)  *FVAF  (1  >1) 

VAi  (J«  2)*FVAR(lB?)*-eVFR(152) 

CALL  VGKAPH(TY^E,VlR,lr  T*  »P»-2»1^»  A  ) 

CALL  VGRA  PH  (TY  J(  I,  FV  17,15*.’, 9,2,1?,  A) 

SR  PETUFN 


.  . 


O'5’.  =  2 


CTN  4.7+476 


1./ 


S»U  "ROUTINE  Gi-4PH(!,TI^C jNKlJN) 

real  kalcorp. 

cor  ‘1C  U  X*l  (152)  » W.dF?)  »FVAV  (1*>2I  ,X1FVAR(t'  2)  ,  XMKVAft(lE2) 
CONMCN/^MkiX/YA.rrMCd  4,2:,  l5j),FTATS(4,3»im  t  ERROS<2,2fl,J:  ) » 

*  TfrbC  <l'c),TY  1  £  Cl-:  ?)  ,  Ve-T  <4,«- ,21> 
nXK£AS10M  A0(1*)  ,-^(l-»)  ,C0d7)  ,13(1M  ,F0(1D  ,£0(17), G0<17>  ,HD(17> 

♦,ocnv  ),:.Dd.') ,30(r),iO(i>) 

DATA  AC  (1)/2.HM“£’<.  EPnr*  +-1  SIGMA/ 

DATA  ;C (31/2. 4  I  viOFT  DYNAMICS  / 


DMA.  AC(  )/2.H  X  CHANNEL  / 

0*1  A  A-:(;  )/2>  M  Fl’K  FOV  / 

fVli  ACO/2..H  TJMr  <«,FC)  / 

DA  1  A  AD(li)/2  i H  PIXELS  * 


CATA  ACC3)/??H  X  CHANNEL  TYWAMTCS  ER-sOF  (S/N=  )  / 

an  d)  -oaisHti)  and  j 
5')  (2)  -cJ(2)=r  )(?1  =  *0(?) 

"0  (7)  =  C0  (•/)=!")(■”)  =:r)(^)  =  FD(7)=G'?(-)=KD(  .*)rAD(7) 

CO  ( o)  =LO(6)=JJO)  =::0(«n=F0(3)=G0(f  )  =HQ ( 8 )  =  .*0(6) 

e  J  (5)800(4) sCJdl)  =  E0  (<n*FD  (9)  sGOmsHCHx)  =  0C  (•»)=?  0  (‘3)  =S0(0)*T 0(3)  = 
J- » .)  (d 

co (*:  i*cn*i-»)<m«t '•)  =Fndu)=FD(i:)=;r»d‘i)sHr(i  )  =oo(1c)=ro<i.  )  = 

1  SJ  d'  >=TJ<1-)=AT(1.) 

*3  (u )  =c> (i;)=''"|di)=r''di)=FD(ii>=;r  di>=*D<u  )*oo(ii)sfn(ii>  = 

'90  (i:)sT  j(n)s;n<:i  > 

FJ  (12  )  -Cu  (l2)  =00  (  12)  =£C(l2)sF3(l?)  =  ;C  (12)eMn(i2)sOD(12)s?.D(12)  = 

-SJ  (1?)=TC  (l2)aA0d.?) 

83  (-  )  rtTjCi)  =  G0(’?>  =  M  (r)  =  RUd)=n<r> 

90  (t  )  =  fcD(d=G0(  >)  =03  (A)  =»/0(l>)  s40(  f  ) 

DMA  CC(i)/2iH  Y  CAIINSL  / 

03  l.)=FD(-)  =  Hj{  I)  =SO(r)=TO(i  )  =C0(  d 
00  (d=FD(f  >  =  H3(A>  =  S)  (d  =  TO(fc)=C*M  f ) 

0A1A  6C(3)/2CH  vriOCITY  / 

D!MJ)*hO(7)sTOO  =*’')(») 
m  {■  )  =»*?(<»)=T3d)  =8*1(4) 

CO  (2)  ■  C!1<7)*S3(*>  =*0(7) 

CO  (~)  =C0(«)=iJ(3)  =  MU ) 

DMA  8C(l3)/‘  1H  X  CHANNEL  /ELJCITY  Eft ’J*  (S/N*  )  / 

UTA  CPd)/2H  ACTUAL  VS.  FILTER  / 

0*1A  CC(2)/2  H  S’GMA  / 

0(1)  =  SC(1)  =  T  0(1)  =  00(1) 

FO  (2)  =  S0<2)  *  r-»C  ?)  =  00(2) 

ro(3)  =  so(3)  =  :nm  =  oo(3) 

O  (*  >  =  so f'»)  =  ;n(  »)  =  Q0(»: ) 

rAT A  CC(13)/e)H  Y  CHAuNrL  TYNAM12E  FR-Oft  (S/N=  )  / 

Da*. A  CCd2)/‘)H  Y  CHANNEL  VELOCITY  ER\Oft  (S/N=  )  / 

OA  *»  A  C0(1»)/‘*|M  FTllt.<  VS.  ACTUAL  SIGid  PLOT  (S/N  =  .  )/ 

DATA  LCd)/2:  ^“^'  T  RACKING  EK\02  / 

DATA  EC(3>/2  M  rCC  dllER  MOOEL  / 
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Performance  Piets  for  Cases  1-12 

This  appendix  contains  the  plotted  outputs  for  cases 
1-12  which  studied  the  performance  of  the  previous  filter 
when  subjected  to  input  parameter  mismatches.  For  most 
cases  four  plots  are  included: 

-  mean  error  ±1  sigma  for  the  filter  estimate  of 
atmospheric  jitter 

-  mean  error  ±1  sigma  for  the  filter  estimate  of  tar¬ 
get  position 

-  real  and  filter-indicated  standard  deviation  of 
atmospheric  jitter 

-  real  and  filter-indicated  standard  deviation  of 
target  position 

Because  of  the  similarity  in  results  between  the  horizontal 
(x)  and  vertical  (y)  channels,  only  the  plots  for  the  hori¬ 
zontal  channel  are  presented. 


Ot'O  03 'O  OO'O  03*0- 

S13XId 

X  CHANNEL  ATMOSPHERICS  ERROR 

Figure  lc.  Case  1  Performance  Plot 


5.60 


Ot'O  06*0  02'0  01  *0 

S~l3XId 


FILTER  V5.  ACTUAL  3IW1R  PLOT 

Figure  If.  Case  1  Performance  Plot 
91 


nrnospHERic  dynamics 
X  CHANNEL 


VWOIS 

SIGMA  CONVERGENCE 


Figure  2a.  Caae  2  Performance  Plot 
92 


4.00  7.00  10.00  13.00  16.00  19.00  22.00 

RUN 


fr-* 

»  m 

o 

o 

09  ;0  00  :0 

09*6- 

00*T-° 

S13XId 

00*0  03*0-  00  *  T-° 

S"13XId 


JRMICS  ERROR 
Performance  Plot 


Figure  2e.  Case  2  Performance  Plot 
96 


ACTUAL  VS.  FILTER 
SIGMA 


SIGMA 

TARGET  DYNAMICS 
Y  CHANNEL 
FLIR  FOV 


09*0  09*0  0**0  02*0  00*0 

VWSIS 


SIGMA  CONVERGENCE 

Figure  2h.  Case  2  Performance  Plot 

99 


4.00  7.00  10.00  13.00  16.00  19.00  22.00 

RUN 


Y  CHANNEL  ATMOSPHERICS  ERROR  CS/N  =  10) 


Figure  2i.  Case  2  Performance  Plot 
100 


19 ‘S 


06'0 


Of'  0 


01*0- 

S13XId 


09*0- 


•  n 

01  *T-C 


r 

Figure 


CHANNEL  DYNAMICS  ERROR 

2j.  Case  2  Performance  Plot 
101 


FILTER  VS.  RCTURL  S1W1R  PLOT  CS/N 


Figure  2k.  Case  2  Performance  Plot 

102 


08  '0 

09'0  0**0 

02  *0 

°o 

mm 

o 

o 

S13XI d 

TURL  SIGWf  pL0T  (S/N 


Performance  Plot 


MEAN  ERROR  «~1  SI  «1  A 
ATHOSTHERIC  OYNAMCS 
X  CHANNEL 
FLIR  FOV 


ft 


0.00  0.80  1.60  2.40  3.20  4.00  4.80  5.60 

TIME  (SEC) 


JM  ERROR  «~i  sicnn 

ruRccr  omomcs 


r 


I 


5.61 


flCTUOL  VS.  FILTER 
SIGMA 


8 


FILTER  VS.  RCTUflL  SIWlfl  PLOT  (S/N  =  101 


Figure  3c.  Case  3  Performance  Plot 
106 


S.6I 


5.61 


00*1 


09*0 


00*0 

S13XId 


09*0- 


00*T-° 


X  CHANNEL  OYI 
Figure  4b.  Case 


4  Performance  Plot 
109 


FILTER 


2  2  1° 
U  U  V- 

a!  x  8 

2  S 


08'0 

09'0  OT’O 

03  ‘0 

00  *0°  1 

S13XId 

FILTER  VS.  ACTUAL  SIGMA  PLOT  ( S/N  =  10) 
Figure  4d.  Case  4  Performance  Plot 

111 


9  *S 


5.6i 


0.00  0.80  1.60  2.40  3.20  4.00  4.80  5.60 

TIME  (SEC) 


riLroi 


FILTER  VS.  qCTdflL  5I0MQ  PLOT  (S/N  =  1  ) 
Figure  5d.  Case  5  Performance  Plot 

115 


0.00  0.80  1.60  2.40  3.20  4.00  4.80  5.60 

TIME  (SEC) 


MEAN  ERROR  *~l  SIGMA 
ATMOSTHERIC  DYNAMICS 


5.61 


now  ckkok  *~i  stem 

TflUGCT  DYNAMICS 


X  CHANNEL  DYNAMICS  ERROR 
Figure  6b.  Case  6  Performance  Plot 
117 


9*S 


HCrUflL  VS.  FILTER 
sicnn 


6c.  Case  6  Performance  Plot 
118 


5.61 


09  ‘  l 


02  *  t 


09  '0 

ST3XI  d 


Ot'O 


00*0° 


FILTER  VS.  qCTJflL  SIGM9  PLCT 
Figure  6d.  Case  6  Performance  Plot 

119 


r 


5.61 


FILTER 


a  z  z 

12  §  I  ° 

w  o  *- 

d  x  a 
2  S 
K 


00*8 

00*9  00** 

00*2 

00*0° 

ST3XI d 

FILTER  VS.  ACTUAL  SIGMA  PLOT  (S/N  =  10) 
Figure  7d.  Case  7  Performance  Plot 

123 


5.6i 


4-1  Slow 


5.61 


IC  DTNRMICS 


09 'I  0  2*1  08*0  0t*0  00 

S13XId 

FILTER  V8.  RCTURL  5IOHR  PLOT  IS/N  =  10) 

Figure  8c.  Case  8  Performance  Plot 

126 


^.00  0.80  1.60  2.40  3.20  4.00  4.80  5.60 

' _  TIME  (SEC) 


09*1  02*t  08*0  oro 

snaxid 

FILTER  V9.  ACTUAL  S1CMA  PLOT  (S/N 

Figure  8d.  Case  8  Performance  Plot 


5.61 


t»W  ERROR  ♦— 1  SIOIR 


Y  CHANNEL  DYNAMICS  ERROR  (8/Ns  1) 

Figure  9b.  Case  9  Performance  Plot 

129 


19  *S 


5.6< 


00  :t 

03  !0  00  !0 

09*6- 

00 *T-° 

siaxid 

MEAN  ERROR  *~1  SIGMR 
TBRCET  OrNRMICS 


X  CHANNEL  DYNAMICS  ERROR 


Figure  10b.  Case  10  Performance  Plot 


133 


5.61 


FILTER  VS.  ACTUAL  S10HA  PLOT 
Figure  10c.  Case  10  Performance  Plot 

13^ 


^J.OO  0.80  1.60  2.40  3.20  4.00  4.80  5.60 

• _  TIME  (SEC) 


/ 


FILTER  VS.  ACTUAL  SIGMA  PLOT 
Figure  lOd.  Case  10  Performance  Plot 


135 


^.00  0.80  1.60  2.40  3.20  4.00  4.80  5.60 

TIME  (SEC) 


HCflN  CltltOR  4-1  SIGMA 

atmospheric  ommics 


X  C 
e  1 


0SPHER1CS  El 
.1  Performarv 
36 


I 


19  S 


00*0 

S13XId 


0**0- 


NflMICS  ERROR  C5/N=  10) 
11  Performance  Plot 


08*0 

09’0  ot*o 

03*0 

00  ’tf3  • 

S13XId 

TRRGET  ommics 


:turl  sigma  plot  (5/n  =  ioj 

LI  Performance  Plot 
L  39 


5.61 


00‘t 

00*3  00*0 

00*2- 

00 • 

S13XId 

MEAN  ERROR  *~l  SIGMR 


X  CHANNEL  0YNAH1CS  ERROR  CS/N=  10) 


19  *S 


5.6i 


OO’t 


oo'e 


oo  *i 


00*2 

S13Xld 

FILTER  V5.  ACTUAL  5I0MA  PLOT  (S/N 
Figure  12d.  Case  12  Performance  Plot 

143 


Appendix  J 


Performance  Plots  for  Cases  13-31 


This  appendix  contains  the  plotted  outputs  for  cases 
13-31  which  studied  the  performance  of  the  previous  filter 
when  subjected  to  the  new  truth  model  changes  formulated  in 
Chapter  II.  For  most  cases  four  plots  are  included: 

-  mean  error  ±1  sigma  for  the  filter  estimate  of 
atmospheric  jitter 

-  mean  error  ±1  sigma  for  the  filter  estimate  of  tar¬ 
get  position 

-  real  and  filter-indicated  standard  deviation  of 
atmospheric  jitter 

-  real  and  filter-indicated  standard  deviation  of 
target  position 

Where  similarities  in  results  existed  between  the  horizontal 
(x)  and  vertical  (y)  channels,  only  the  plots  for  the  hori¬ 
zontal  channel  are  presented. 
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Appendix  K 

Input  Parameters  for  Cases  32-37 


Refer  to  List  of  Symbols  for 
•Description  of  Variables 
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Appendix  L 

Performance  Plots  for  Cases  32-37 

This  appendix  contains  the  plotted  outputs  for  cases 
32-37  which  studied  the  performance  of  the  new  filter  formu¬ 
lated  in  Chapter  IV  and  discussed  in  Chapter  V.  For  most 
cases  four  plots  are  included: 

-  mean  error  ±1  sigma  for  the  filter  estimate  of  tar¬ 
get  position 

-  mean  error  +1  sigma  for  the  filter  estimate  of  tar¬ 
get  velocity 

-  real  and  filter-indicated  standard  deviation  of 
target  position 

-  real  and  filter-indicated  standard  deviation  of 
target  velocity 

Where  similarities  in  results  existed  between  the  horizontal 
(x)  and  vertical  (y)  channels,  only  the  plots  of  the  hori¬ 
zontal  channel  are  presented. 
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